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Abstract 

Current observations of double neutron stars provide us with a wealth of infor- 
mation that we can use to investigate their evolutionary history and the physical 
conditions of neutron star formation. Understanding this history and formation con- 
ditions further allow us to make theoretical predictions for the formation of other 
double compact objects with one or two black hole components and assess the de- 
tectability of such systems by ground-based gravitational-wave interferometers. In 
this paper we summarize our group's body of work in the past few years and we 
place our conclusions and current understanding in the framework of other work in 
this area of astrophysical research. 
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1 Introduction 



Gravitational waves were theoretically predicted many decades ago by Ein- 
stein and his theory of General Relativity, and their existence has been in- 
directly _confixmed__bxthe_pioneering observations of the Hulse- Taylor binary 
pulsar (IHulse fc Taylorll 19751 ). In the coming years, the physics and astronomy 
community anticipates the first direct detection of gravitational waves (GW) 
and the birth of a new research field: GW astrophysics, where observations 
with the first ground-based interferometers (such as LIGO, GEO, and Virgo) 
contribute to the astrophysical understanding of the sources. Double neutron 
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stars (DNS) and their inspiral due to GW emission represent the primary tar- 
get source. This is primarily because they are kno wn to exist in the G alaxy, 
with the Hulse- Taylor binary and the binary pulsar (IBurgay et al.ll2003l ) being 
the staples of this class of binary compact objects. Double compact objects 
with black hole (BH) components (i.e., BH-NS and BH-BH) have not yet 
been discovered in nature, but our current understanding of binary star evo- 



ble time ( 


Belczvnski et al. 


2002b|: 


Brown 


1995; 


Porteeies Zwart & Yuneelson 


1998; 


O'Shauehnessv et al. 


2005b; 


Hurlev et al.ll2002l: 


Porteeies Zwart & Verbunt 


1996; 


Frver et al. 


1998; 


Bethe & Brown 1999h. 



The sample of known DNS in the Galaxy has been very well studied and 
monitored since their discovery and as a result there is a wealth of information 
about their physical properties. This observational sample provides us with 
a unique opportunity to investigate the evolutionary history of these systems 
and a number of studies have addressed this basic question from different 
points of view. Another way to address the question of the formation of double 
compact objects is to develop population synthesis models that follow the 
evolution of large ensembles of binaries. Such models allow the detailed study 
of evolutionary channels that lead to double compact object formation and 
assess their relative efficacy and the main binary evolution phases that affect 
the final properties of the binaries. 



The calculation of inspiral rates for double compact objects and associated 
event rates for ground-based GW detectors lies at the center of research ac- 
tivity rela ted to GW sear ches in the last decade. Th ese calculations were pio- 
neered by iPhinneyl ( Il99ll ) and lNarayan et al.l (Il99ll ). Since then, a number of 
studies have considered this problem and revised these initia l resul ts, following 
the methods introduced in 1991. A few years ago lKim et al.l (120031 ) introduced 
a novel statistical method for the calculation of inspiral rates based on the ob- 
served systems (empirical estimates); this method allowed us for the first time 
to associate rate estimates with statistical and systematic uncertainties and 
to place these estimates on the basis of confidence levels. Inspiral rates have 
also been computed based on population synthesis cal culations that explore 



a limite d number of model assumptions. More recently lO'Shaughnessy et al 



( I2005bl ) have combined the two methods and have used the empirical DNS 
estimates as constraints on model predictions, especially for BH binaries on 
which we have no direct guidance from observations. 



In this paper we summarize our group's work in this area in the context of 
double compact object formation, their properties and inspiral rates and in 
connection to other studies that have appeared in the literature. 
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2 The Recent Evolutionary History of Known Close Double Neu- 
tron Stars 



According to our current knowledge (for more details see 


Bhattacharva k van den Heuvel 


1991 




Brown 


1995 


; Belczvnski et al. 


2002b: 


Tauris k van den Heuvel 


2006 





double neutron stars (DNSs) form from primordial binaries in which, possibly 
after some initial mass-transfer phase, both component stars have masses in 
excess of ~ 8 — 12 M Q . After the primary explodes in a supernova (SN) explo- 
sion to form the first neutron star (NS), the binary enters a high-mass X-ray 
binary phase in which the NS accretes matter from the wind of its companion. 
The phase ends when the companion star evolves of the main sequence and 
engulfs the NS in its expanding envelope. The NS then spirals in towards the 
core of the companion until enough orbital energy is transferred to the en- 
velope to expel it from the system. When the envelope is ejected, the binary 
consists of the NS and the stripped-down helium core of its former giant com- 
panion, orbiting each other in a tight orbit. If the NS is able to accrete during 
its rapid inspiral, a first "recycling" may take place during which it is spun up 
to millisecond periods. The helium star then evolves further until it, in turn, 
explodes and forms a NS. Depending on the mass of the helium star and the 
size of the orbit at the time of the explosion, the SN may be preceded by a 
second recycling phase when the helium star fills its Roche lobe and tr ansfers 
mass to the NS flBelczvhski k Kalogerall200ll ; iBelczvnski et al.ll2002al lbh. 



In what follows in this section we use this basic evolutionary framework and 
the complete, up-to-date set of observational constraints (at the time of this 
writing) on three DNS Galactic systems to investigate the physical properties 
of progenitor binaries at the time of the second SN exp losion. We primarily 



summ arize the work presented in IWillems et al.l (120061 ) and IWillems et al. 



( 120041 ) with many comparisons to other published studies addressing similar 
questions. 



2.1 PSR J 01 31 -3039: The Double Pulsar 



The recent discovery of the strongly relativistic binary pulsar (IBurgay et al 



2003 ) which is also the first eclipsing double pulsar system found in our Galaxy 



( ILyne et al.ll2004l ) has res parked the interest in the evolutionary history and 
form a tion of DNS systems (IWillems k Kaloger al l2004J : iDewi k van den Heuvel! 



2004: IWillems et all 1200411 S hortl y after the discovery of P SR J0737-3039, 



Dewi k van den Heuvel! (120041 ) and lWillems k Kalogeral (120041 ) independently 



derived that, right before the second SN explosion, the binary was so tight that 
the helium star must have been overflowing its critical Roche lobe. This con- 
clusion, for the first time, strongly confirmed the above formation channel 
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for DNS binaries. The observed 22.7ms pulsar (hereafter PSR J0737-3039A or 
pulsar A) then corresponds to the first-born NS, and its 2.8 s pulsar companion 
(hereafter PSR J0737-3039B or pulsar B) to the second-born NS. 

The main properties of the double pulsar system and it s component stars re l- 
evant to this investigation are summarized in Table 1 of I Willems et al.l (120061 ). 



The kine matic properties of PS R J0737-3039 have undergone significant revi- 



sion since iRansom et al.l (120041 ) used interstellar scintillation measurements 
and inferre d a velocity comp onent perpendicular to the line-of-sight of ~ 
141kms _1 . I Coles et al.l (120051 ) derived a reduced velocity of ~ 66kms _1 by 
incorporating anisotropies of the interstellar medium in the interstellar scintil- 
lation model. This reduced veloc i ty, ho wever, strongly depends on the adopted 
anisotropy model. iKramer et al.l (120051 ). on the other hand, used pulsar timing 



measurements to derive a firm model-independent upper limit of 30 km s 1 on 
the transverse velocity. 

The kick imparted to pulsar B at birth is expected to tilt the orbital plane 
and misalign pulsar A's spin axis wit h respect to the post-SN orbital angular 
momentum axis (e.g.. lKalogerall2000l ). The degree of misalignment depends on 
both the magnitude and the direction of the kick, and therefor e yields a valu- 
able p iece of information on pulsar B's natal kick velocity. In I Willems et al. 
(120041 ) we showed that for the pre- and post-SN orbital parameters compati- 
ble with all available observational constraints, misalignment angl e s betw een 
approximately 70° and 110° are highly unlikely. Manchester et al.l (120051 ) de- 
rived observational constraints on the spin-orbit misalignment based on the 
stability of pulsar A's mean pulse profile over a time span of 3 years. The 
authors concluded the angle to be smaller than ~ 60° or larger than ~ 120°, 
in agreement with the theoretic al predictions of IWillems et al.l (120041 ) . These 
misalignment limits inferred by Manchester et al.l (120051 ) are incorporated in 
the list of available observational constraints. 



Among the constraints inferred from observations, the most uncertain param- 
eter affecting the reconstruction of the system's formation and evolutionary 
history is its age. Characteristic ages, defined as half the ratio of the spin 
period to the spin-down rate, are the most c ommonly adopted ag e estimators, 
but are known to be quite unreliable (e.g. IKramer et al.l 120031 ) . In the case 
of PSR J0737-3039, the spin evolution of pulsar B is furthermore ve ry likely 
affected by torques exerted b y pulsar A on pulsar B (jLyutikovl 120041 ). adding 
to the uncertainties of its age. lLorimer et al.l ( 20051 ) therefore derived an alter- 
native age estimate by noting that, according to the standard DNS formation 
channel, the time expired since the end of pulsar A's spin-up phase should 
be equal to the time expired since the birth of pulsar B. By combining this 
property with a selection of different spin-down models, the authors derived a 
most likely age of 30-70 Myr, but were unable to firmly exclude younger and 
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older ages. A third age estimate can be obtained by assuming that pulsar A 
was recycled to the maximum spin rate and calculating the time required for 
it to spin down to the currently observed value. This so-called s pin- down age 



provi des an upper limit to the age of the system of 100 Myr (IBurgay et al. 



20031 ). In view of these uncertainties, we derive constraints on the formation 



and evolution of PSR J0737-3039 for three different sets of age assumptions: 
(i) < t < 100 Myr, (ii) 30 < r < 70 Myr (the most likely age range de- 
rived by lLorimer et al.l (120051 )). and (iii) r ~ 50 Myr (the characteristic age of 
pulsar B). 



2.1.1 Basic assumptions and outline of the calculation 

In the following subsections we use the presently known observational con- 
straints to reconstruct the evolutionary history of PSR J0737-3039 and con- 
strain its properties at the formation time of pulsar B. More specifically, our 
aim is to derive a probability distribution function (PDF) for the magnitude 
and direction of the kick velocity imparted to pulsar B at birth, the mass of 
pulsar B's progenitor immediately before it explodes in a SN explosion, and 
the orbital separation of the component stars ri ght before the SN explo sion. 



Other recen t stud ies with similar goals include iPiran fc Shavivl (120061 ) and 



Stairs et al. ( 20061 ). Adopting the standard DNS formation channel, the pro- 



genitor of the double pulsar right before the formation of pulsar B, consists 
of the first-born NS, pulsar A, and the helium star progenitor of the second- 
born NS. Since the formation of the second NS is preceded by one or more 
mass-transfer phases, tidal forces can safely be assumed to circularize the or- 
bit prior to the formation of pulsar B. In what follows, we refer to the times 
right before and right after the formation of pulsar B by pre-SN and post-SN, 
respectively. The analysis presented consists of four major parts. 

Firstly, the motion of the system in the Galaxy is traced back in time up 
to a maximum age of 100 Myr. The goal of this calculation is to derive the 
position and center-of-mass velocity of the binary right after the formation 
of pulsar B, and use this information to constrain the kick imparted to it at 
birth. In order to determine possible birth sites, our analysis is supplemented 
with the assumptions that the primordial DNS progenitor was born close to 
the Galactic plane, and that the first SN explosion did not kick the binary 
too far out of it. The first assumption is conform with our current knowledge 
that massive stars form and live close to the Galactic plane (their typical scale 
height is ~ 50-70 pc). The second assumption on the othe r hand neglects a 



small fraction of system s formed with high space velocities (jPfahl et al.l 12002 



Belczynski et al.ll2002bl ). Under these assumptions, the binary is still close to 
the Galactic plane at the formation time of pulsar B. Possible birth sites can 
thus be obtained by calculating the motion in the Galaxy backwards in time 
and looking for crossings of the orbit with the Galactic mid-plane. The times 
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in the past at which the plane crossings occur yield kinematic estimates for 
the age of the DNS, while comparison of the system's total center-of-mass 
velocity with the local Galactic rotational velocity at the birth sites yields the 
system's post-SN peculiar velocity. 



Secondly, the orbital semi-major axis and eccentricity right after the SN explo- 
sion forming pulsar B are determined by integrating the equations governing 
the evolution of the orbit due to gravitational wave emission backwards in 
time. The integration is performed for each Galactic plane crossing found 
from the Galactic motion calculations. Since each crossing yields a different 
kinematic age, and thus a different endpoint of the reverse orbital evolution 
calculation, the post-SN parameters are a function of the considered Galactic 
plane crossing. 



Thirdly, the conservation laws of orbital energy and angular momentum are 
used to map the post-SN binary parameters to the pre-SN ones. The map- 
ping depends on the kick velocity imparted to pulsar B at birth and results in 
constraints on the magnitude and direction of the kick velocity, the mass of 
pulsar B's pre-SN helium star progenitor, and the pre-SN orbital separation. 
The pre-SN orbital eccentricity is assumed to be zero, as expected from strong 
tidal forces operating on the binary during the mass-transfer phase (s) responsi- 
ble for spinning up pulsar A. The kick and pre- and post-SN binary parameters 
are then furthermore subjected to the requirements that they be consistent 
with the post-SN peculiar velocity obtained from the Galactic motion calcu- 
lations and with the observationally inferred spin-orbit misalignment angle of 
pulsar A. The latter constraint requires the additional assumption that tidal 
forces align the pre-SN rotational angular momentum vector of pulsar A with 
the pre-SN orbital angular momentum vector. 



Fourthly, the solutions of the conservation laws of orbital energy and angu- 
lar momentum are weighted according to the likelihood that they lead to the 
system's currently observed position and velocity in the Galaxy. A PDF of 
the admissible kick velocity and progenitor parameters is then constructed by 
binning the solutions in a multi-dimensional "rectangular" grid. The maxi- 
mum of the resulting PDF yields the most likely magnitude and direction of 
the kick velocity imparted to pulsar B at birth, mass of pulsar B's pre-SN 
helium star progenitor, and pre-SN orbital separation. We also investigate the 
sensitivity of the PDF to the adopted assumptions by systematically varying 
the underlying assumptions, such as, e.g., the age and the magnitude of the 
transverse systemic velocity component. 
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2.1.2 Kinematic age and history 



Tracing the motion of PSR J0737-3039 back in time requires the knowledge of 
both its present-day position and 3-dimensional space velocity. Unfortunately, 
no method is presently available to measure radial velocities of DNSs, so that 
the knowledge of their space velocity is limited to the component perpendicular 
to the line-of-sight (transverse velocity). At the time of this analysis, only an 
upper limit was available for the double pulsaJ"*"! The direction of the motion 
in the plane perpendicular to the line-of-sight is still unknown. Therefore we 
explore the system's kinematic history in terms of two unknown parameters: 
(i) the radial component V r of the present-day systemic velocity, and (ii) the 
orientation Q of the transverse velocity in the plane perpendicular to the line- 
of-sight. Since the results presented i n this paper do not d epend on the exact 



definition of Q, we refer the reader to lWillems et al.l (120041 ) for a more detailed 



discussion and definition. We calculate the motion of the system in the Galaxy 



backwards in time using the Galactic poten tial of ICarlberg Sz Innanenl (119871 ) 



with updated model parameters derived by iKuijken fc Gilmord (119891 ). Since 
at the time of the analysis kinematical constraints provided only an upper 
limit of SOkms" 1 on the transverse systemic velocity, and the calculation of 
the past orbit requires precise starting values for both the present-day position 
and the present-day velocity, we calculate the motion backwards in time for 
two specific velocities: V t = lOkms -1 and V t = 30kms _1 in what follows. 

The number of Galactic plane crossings associated with each V r and Q, within 
the adopted age limit of 100 Myr, can be anywhere between 1 and 5. The times 
in the past at which the system crosses the Galactic plane furthermore yield 
kinematic estimates for the age of the DNS, while subtraction of the Galactic 
rotational velocity from the total systemic velocity at the birth sites yields the 
peculiar velocity right after the formation of pulsar B. 

Since there are no a priori constraints on the magnitude of the radial velocity, 
we examine the effect of incorporating large radial velocities in the analysis by 
weighting each considered radial velocity according to a pre-determined radial 
velocity distribution. For this purpose, we performed a population synthesis 
study of Galactic DNSs, inc luding their kinematic evolution in the potential of 



Carlberg fc Innanenl (119871 ). Theoretical radial velocity distributions are then 
obtained by taking a snapshot of the population at the current epoch and de- 
termining the radial velocity for each DNS in the sample. The resulting PDFs 
are found to be represented well by Gaussian distributions with means of 
Okms" 1 and velocity dispersions of 60-200 km s -1 , depending on the adopted 
population synthesis assumptions. For comparison, we also consider a uniform 
radial velocity distribution in which all radial velocities between -1500kms _1 
and 1500kms _1 are equally probable. For the other unknown parameter, the 

1 At present a measurement of 10 km s _1 has been reported by Kramer et al. ( 20061 ). 
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Fig. 1. Distribution of post-SN peculiar velocities for a present-day transverse ve- 
locity of 30kms _1 , kinematic ages ranges of 0-100 Myr and 49-51 Myr, and radial 
velocity distributions varying from a uniform distribution (solid line) to Gaussian 
distributions with velocity dispersions of 60kms _1 (long-dashed line), 130 km s -1 
(short-dashed line), and 200 kms -1 (dash-dotted line). For clarity, the PDFs are 
offset from each other by an arbitrary amount. 

orientation angle Q of the transverse velocity component in the plane perpen- 
dicular to the line-of-sight, we assume a uniform distribution between 0° and 
360°. 



In Fig. 1, the distribution of post-SN peculiar velocities obtained by trac- 
ing the motion of PSR J0737-3039 in the Galaxy back in time is shown for 
a present-day transverse velocity component of 30kms _1 , and kinematic age 
ranges of 0-100 Myr and 49-51 Myr. These probability distributions are cal- 
culated considering weights according to (i) the probability that the system 
has a present-day radial velocity V r (assumed to be distributed according to a 
uniform or a Gaussian distribution), (ii) the transverse velocity has a direction 
angle Q (assumed to be uniformly distributed), and (iii) the system is found at 
its current position in the Galaxy (determined by the time the system spends 
near its current position divided by its kinematic age). 

Here, we determine the probability of finding the system at its current po- 
sition in the Galaxy in three dimensions (i.e., besides the vertical distance 
to the Galactic plane, we also consider the radial and azimuthal position in 
the plane). For this purpose, we determine the time the system spends in a 
sphere with radius R centered on its current location for all V r - and ^-values 
considered for the derivation of the possible birth sites. The probability to 
find the system near its current position is then determined as the time it 
spends in this sphere divided by its age. For the latter, we adopt the kine- 
matic ages obtained by tracing the motion of the system backwards in time, 
so that the probability of finding the system near its current position depends 
on the radial velocity V r , the proper motion direction Q, and the Galactic 
plane crossing considered along the orbit associated with V r and Q. We per- 
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formed some test calculations to successfully verify that the results based on 
the thus determined probabilities are insensitive to the adopted value of R, as 
long as it is sufficiently small for the sphere to represent a local neighborhood 
near the system's current position. For the results presented in this paper we, 
somewhat arbitrarily, choose R = 50 pc. 

Ultimately, the probability that, for a given pair of V r and Q, the post-SN 
peculiar velocity of the binary is equal to Vcm is given by 

p(vbM|v r ,n)<x £ \i; > \i(v r ,n), (i) 

where N(V r ,Q) is the number of Galactic plane crossings associated with V r 
and Q, T(V r , Q) is the time the system spends near its current position for 
the orbit associated with V r and Q, and Tki n ,i (V r , Q) is the kinematic age 
corresponding to the i-th plane crossing associated with V r and Q. The factor 
\i (V r , Q) is equal to 1 when the i-th plane crossing along the orbit associated 
with V r and f2 yields a post-SN peculiar velocity equal to Vcm, and equal to 
otherwise. The total probability that the post-SN peculiar velocity of the 
binary is equal to Vcm is then obtained by integrating Eq. (1) over all possible 
values of V r and Q: 

p(v cm )oc J J p(v CM \v r ,n)p(v r )p(n)dv r dn. (2) 

n v r 

Here P(V r ) is the probability distribution of the unknown radial velocity V r , 
and P(Q) the probability distribution of the unknown proper motion direction 
fl in the plane perpendicular to the line-of-sight. 

For ages of 0-100 Myr, post-SN peculiar velocities up to 100 km s -1 are likely 
for all four radial velocity distributions. The highest post-SN peculiar veloc- 
ities are found for the Gaussian radial velocity distributions with velocity 
dispersions of 130 km s^ 1 and 200 km s -1 (Vcm values remain likely up to 200- 
300 km s -1 ) and the uniform radial velocity distribution (Vcm values follow an 
almost flat distribution from 300 km s" 1 to 1500 km s _1 ). For ages of 49-51 Myr, 
the post-SN peculiar velocity distributions all have a peak at 50kms _1 . In the 
case of the uniform radial velocity distribution and the Gaussian distributions 
with velocity dispersions of 130 km s" 1 or 200 km s _1 , additional peaks of al- 
most equal height occur at even higher post-SN peculiar velocities. Similar 
conclusions apply to the post-SN peculiar velocity distributions obtained for 
a present-day transverse velocity component of lOkrns -1 . Despite the small 
lower limit on the present-day transverse systemic velocity, a wide range of 
non-negligible post-SN peculiar velocities is thus possible. We note however 
that the distributions shown here incorporate all possible post-SN peculiar 
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velocities obtained from tracing the motion of the system backwards in time 
and that some of these may require SN kicks and mass loss that are incompat- 
ible with the observational constraints on the orbital elements and component 
masses of the double pulsar. 



2.1.3 Progenitor constraints 



After the formation of pulsar B, the evolution of the system is expected to 
be driven exclusively by the e mission of gravita t ional waves. We use the dif- 
ferential equations derived by lJunker fc Schaferl (119921 ) which are valid up to 
3.5 post-Newtonian order of approximation, and find the resulting post-SN 
orbital separation A to be always between 1.26 R Q and 1.54 R Q , and the post- 
SN orbital eccentricity e between 0.088 and 0.12. Since the kinematic ages 
are functions of the radial velocity V r , the proper motion direction Q, and 
the Galactic plane crossing considered along the past orbit associated with V r 
and Q, the particular values of A and e associated with a given Tki n are also 
functions of these variables. 



The pre- and post-SN binary parameters and the kick velocity imparted to 
pulsar B at birth are related by the conservation laws of orbital energy and 
angular momentum. For a circular pre-SN orbit, the relations can be found in 
(Hills 1 9831; iBrandt fc Podsiadlowskilll995l ; lKalogeralll996l ; iFrver fc Kalogera 
19971 : iKalogera & Lorimerll2000h . 



The requirement that these equations permit real solutions for the kick mag- 
nitude Vfc, the kick orientation angles 9, (f>, the NS progenitor mass M , and 
the pre-SN orbital semi-major axis A , imposes constraints on all of these 
parameters. For a m athematical fo r mulat ion of these constraints, we refer 
to Eqs. (21)-(27) in IWillems et al.l (120051 ), and references ther ei n ( a more 



compact descr i ption can also be found in IWillems fc Kalogeral (12004 ) and 
Willems et al.l (120041 ) . We here merely recall that the constraints express that: 
(i) the binary components must remain bound after the SN explosion, (ii) the 
pre- and post-SN orbits must pass through the instantaneous position of the 
component stars at the time of the SN explosion, and (iii) there is a lower and 
upper limit on the degree of orbital contraction or expansion associated with 
a given amount of mass loss and a given SN kick. 

Besides the changes in the orbital parameters and the mass of the explod- 
ing star, the SN explosion also imparts a kick velocity to the binary's center 
of mass and tilts the post-SN orbital plane with respect to the pre-SN one. 
The center-of-mass velocity and tilt angle add additional constraints on the 
solutions of the conservation laws of orbital energy and angular momentum. 

The constraints on the progenitor of PSR J0737-3039 resulting from the dy- 
namics of asymmetric SN explosions arise solely from tracing the evolution 
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of the current system properties backwards in time. The pre-SN orbital sep- 
arations and pulsar B progenitor masses found this way are, however, not 
necessarily accessible through the currently known DNS formation channels. 
Further constraints on the progenitor of pulsar B can therefore be obtained 
from stellar and binary evolution calculations. 

Firstly, a lower limit on the mass of pulsar B's pre-SN progenitor is given by 
the requirement that the star must be massive enough to evolve into a NS 
rather than a white dwarf. According to our current understanding of helium 
star evoluti on, the minimum helium star mass required for NS formation is 
2.1-2.8 M Q flHabetslll98i iTauris fc van den Heuvelll2006l ). The actual helium 
star minimum mass is, however, still considerably uncertain due to the poorly 
understood evolution of massive stars and possible interactions with close 
binary companions. 



We impose a lower limit of 2.1 M on the mass of pulsar B's pre-SN helium 
star progenitor. However, in the light of recent suggestions that the progeni- 
tor of pulsar B may have been significantly less massive than the conventional 
lower limit of 2.1 M & , we also explore the possibilit y of progenitor masses 
as low as pulsar B's pre sent-day mass of 1.25 M (IPiran fc Shavivl l2005al ; 
Podsiadlowski et al.l 120051 ). Unless our current understanding of helium star 
evolution is seriously flawed, this scenario implies that the progenitor of pul- 
sar B must have lost a significant amount of mass (at least ~ 0.7 M & ) after it 
had already established a sufficiently massive core to guarantee the occurrence 
of a future SN explosion. We note that this possibility is included in the binary 
population synthesis calculations used to derive the adopted theoretical DNS 
radial-velocity distributions. 



Secondly in IWillems fc Kalogeral (120041 ). we showed that the pre-SN binary 
was so tight (1.2 R e < A < 1.7 R & ) that the helium star progenitor of pul- 
sar B m ust have been overflowing its R oche lobe at the time of its SN explosion 
(see also Dewi fc van den Heuvellliool . An upper limit on the progenitor mass 
is therefore given by the requirement that this mass-transfer phase be dynami- 
cally stable (otherwise the components would have merged and no DNS would 
have formed) . Based on the evolut ionary tracks for NS + helium star binaries 



calculated by llvanova et al.l (120031 ). we adopt an upper limit of 3.5 on the mass 
ratio of the pre-SN binary to separate dynamically stable from dynamical ly 
unstable Roche-lobe overflow (see also lDewi et al.l 120021 ; iDewi fc Polsll2003l ). 



2.1.4 The most likely kick velocity and progenitor properties 



The most likely kick velocity and progenitor parameters obtained from the 
procedure outlined i n the previous sections are summarized in Table IV of 
(IWillems et al.l 120061 ) . The most likely pulsar B kick velocity is smaller than 
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50kms _1 only when the transverse velocity component has a magnitude of 
10 km s _1 , the minimum helium star mass required for NS formation is allowed 
to be as low as 1.25 M Q , and the age of the system is assumed to be between 30 
and 70 Myr. All other models yield most likely kick velocities of 50-180 km s _1 . 
When a minimum pre-SN helium star mass of 2.1 M is imposed, the kicks are 
always strongly favored to be directed opposite to the helium star's pre-SN 
orbital motion (most likely cos# ~ —0.90 ±0.05). When the constraint on the 
minimum pre-SN helium star mass is relaxed, the most likely kick direction 
can shift significantly and can even become perpendicular to the helium star's 
pre-SN orbital velocity. Allowing pre-SN helium star masses down to 1.25 M 
furthermore always leads to most likely progenitor masses of 1.3 — 2.0 M Q , 
except when Vt = lOkrns^ 1 , the age of the system is between and 100 Myr, 
and the radial velocities are distributed uniformly or according to a Gaussian 
with a velocity dispersion of 200 km s^ 1 . In the latter cases, the most likely 
progenitor masses are 2.5-2.7 M . 

In summary, while small kick velocities of just a few tens of kms -1 could be 
favored for some models, the majority of the models yields most likely val- 
ues of 50-180 km s -1 . Progenitor masses below 2.1 M Q are furthermore not 
required to explain the system properties, although they are generally favored 
when helium stars below 2.1 M are still assumed to be viable NS progeni- 
tors. We note though that the results presented are based on the assumption 
that all kick velocity magnitudes 14 are equally probable. This assumption is 
inconsistent for models using the Gaussian radial velocity distributions based 
on population synthesis calculations of coalescing DNSs. In particular, these 
calculations all adopt a Maxwellian rather than a uniform kick velocity distri- 
bution. In all cases, weighing the kick velocities according to the Maxwellian 
underlying the derivation of the radial velocity distributions would, however, 
shift the most likely kick velocities to higher V k values. This reinforces our 
conclusion that the presently known observational constraints not necessarily 
disfavor kick velocity magnitudes of 100 km s -1 or more. 



2.1.5 Confidence limits on the kick velocity and progenitor parameters 

For illustration, some representative 1-D PDFs used for the calculation of the 
confidence limits in the case of a present-day transverse velocity component of 
lOkms" 1 are shown in Fig. 2. Panels (a)-(c) show the kick velocity distribu- 
tions resulting from different present-day radial velocity distributions for ages 
ranges of 0-100 Myr and 49-51 Myr, and minimum pre-SN helium star masses 
of 1.25 M Q and 2.1 M Q . For a given age range and minimum pre-SN helium 
star mass, the PDFs show a peak which is most pronounced when a Gaussian 
radial velocity distribution with a velocity dispersion of 60kms _1 is consid- 
ered, and which widens with increasing radial velocity dispersions. When the 
age is assumed to be 0-100 Myr and the minimum pre-SN helium star mass 
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Fig. 2. One-dimensional PDFs illustrating some of the dependencies of the derived 
pulsar B kick velocity and progenitor properties on the adopted model assumptions. 
All plots are for a present-day transverse velocity of 10 km s _1 . Solid lines correspond 
to uniform radial velocity distributions and dashed lines to Gaussian distributions 
with velocity dispersions of 60kms -1 (long-dashed lines), 130 km s -1 (short-dashed 
line), and 200kms _1 (dash-dotted line). For clarity, the PDFs are offset from each 
other by an arbitrary amount. Panels (a)-(c) show the distributions of kick velocity 
magnitudes panels (d)-(e) the distributions of kick direction cosines cos#, and 
panel (f) the distributions of pre-SN helium star masses Mq. 



13 



1.25 Mq, there is a clear tendency for the 1-D PDFs to favor kick velocities 
of 50kms _1 or less (although this becomes significantly less pronounced with 
increasing radial velocity dispersions). This trend shifts towards favoring kick 
velocities of 50-100 km s _1 when the age range is narrowed to 49-51 Myr. More- 
over, when the age range is kept fixed at 0-100 Myr, but the minimum pre-SN 
helium star mass is increased to 2.1 M Q , the favored range of kick velocities 
shifts to 100-150 km s _1 . In the latter case, the kick ve locity is furthermore 



always required to be larg er than ~ 60 km s 1 (see also I Willems fc Kalogera 



2004J : IWillems et alJl2004[ ). 



Panel (f), finally, shows the distribution of possible pre-SN progenitor masses 
for different present-day radial velocity distributions, DNS ages of 0-100 Myr, 
and a minimum pre-SN helium star mass of 1.25 M Q . The distributions all 
favor progenitor masses of 1.4-1.5 M , with the preference for this mass range 
being strongest for present-day radial velocity distributions with small radial 
velocity dispersions. Distribution functions for a minimum pre-SN helium star 
mass of 2.1 M look practically the same as the ones displayed panel (f) if 
they were cut off at 2.1 M . 



2.1.6 Summary and concluding remarks 



One of the results of the above analysis is that marginalizing the full five- 
dimensional progenitor PDF to 1-D or 2-D distributions for the pulsar B 
kick velocity and progenitor mass has a major effect on the determination of 
their most likely values. When the full multi-dimensional PDF is examined, 
it is clear that although some sets of prior assumptions indeed favor low kick 
velocities and low progenitor masses the majority of the models favor kick 
velocities of 50-180 km s _1 and progenitor masses of 1.45-2.75 M Q . 



In particular, if the transverse systemic velocity is assumed to be lOkm/s and 
helium stars less massive than 2.1 M are assumed to be viable NS progeni- 
tors, the most likely progenitor mass can vary from 1.35 M & to 2.65 M Q , de- 
pending on the assumed systemic age and radial velocity distribution. Hence, 
whether progenitor masses greater than 2.1 Mq are statistic ally likely or un- 
likely depends strongly on the adopted prior assumptions (cf. IPiran &: Shaviv 
2005al lbl). Most likely progenitors with Mq < 2.1 can furthermor e also be 



associated with kick velocities of up to 100kms~ 1 . IStairs et al.l (120061 ) recently 
ar rived at similar concl usions based on updated proper motion measurements 
by lKramer et al.l (120061 ) . The updated measurements yield a proper motion of 
lOkm/s in a direction nearly parallel to the plane of the Galaxy. 



We also find that the proximity of PSR J0737-3039 to the Galactic plane and 
the small proper motion do not pose stringent constraints on the kick velocity 
and progenitor mass of pulsar B. Instead, the constraints are predominantly 
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determined by the orbital dyn amics of asymm e tric SN ex plosions. This is 
i n cont rast to the conclusions of iPiran fc Shavivl (j2005al ) and iPiran fc Shaviv 
( 12005bl ) who emphasize that the proximity of the double pulsar to the Galactic 
plane implies that pulsar B most likely received only a small kick at birth and 
that its progenitor most likely had mass of ~ 1.45 M & . 



Hence, based on the currently available observational constraints, a wide range 
of progenitor and kick velocity properties are favored for PSR J0737-3039B. 
In particular, the constraints are compatible with a conventional hydrody- 
namical or ne utrino-driven SN explosion from a helium sta r more massive 



than 2.1 M & (IHabets! Il986l ; iTauris fc van den Heuvell 120061 ). as we ll as an 



electron-cap t ure SN from a helium star less massive than 2.5 M & (INomoto 



1984 119871 ). iPodsiadlowski et al.l (120051 ) have speculated that the electron- 



capture SN mechanism may be typical for close binaries and that it may be 
accompanied by much smaller kicks than hydrodynamical or neutrino-driven 
SN explosions. Consequently, if pulsar B is formed through an electron capture 
SN and if electron capture SNe are accompanied by small kicks, the mass of 
pulsar B's pre-SN helium star progenitor must be smaller than 2.1 M & (oth- 
erwise the kick is always larger than ~ 60kms _1 ). Since it is unlikely that 
future observations will lead to new constraints on the smallest possible pul- 
sar B progenitor mass, further insight to the formation mechanism of pulsar B 
should be sought in core-collapse simulations of low-mass (< 2.1 M ) helium 
stars and population synthesis studies of PSR J0737-3039-type binaries and 
their progenitors. 



2.2 PSR B 1534+1 2 



The relativistic binary radio pulsar PSRB1534+12 too has an accurately mea- 
sured proper motion with a known direction in right ascension and declination, 
so that its kinematic history and progenitor constraints depend only on the 
unknown radial velocity V r . In order to derive these constraints, we adopt the 
spin-down age T& = 210 Myr as an upper limit to the age of the system. The 
other physical parameter s of PSR B 1534+ 12 relevant to the derivation are 
summarised in Table 1 of IWillems et al.l (120041 ) . 



Following the same arguments as for PSR J0737-3039, we assume that the 
progenitor of PSRB1534+12 was close to the Galactic plane at the time of 
the second SN explosion and that its pre-SN systemic velocity was almost 
entirely due to Galactic rotation. The possible birth sites of the DNS are then 
obtained by tracing the Galactic motion of the system backwards in time as a 
function of the unknown radial velocity V r . We find that, within the imposed 
age limit of 210 Myr, no Galactic plane crossings occur when V r > 250 km s -1 , 
whereas up to four crossings may occur when V r < 250 km s -1 . Since four 
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Fig. 3. Variations of the kinematic age (left-hand axis) and post-SN peculiar velocity 
(right-hand axis) of PSRB1534+12 as a function of the unknown radial velocity V r , 
for cases A, B, and C. The kinematic ages are represented by red lines and the 
peculiar velocities by blue lines. 

disk crossings only occur for relatively few and rather fine-tuned Galactic 
trajectories, we leave these cases aside and focus on the possible birth sites 
associated with the first (case A), second (case B), and third (case C) Galactic 
plane crossings. 

The kinematic ages and post-SN peculiar velocities associated with the Galac- 
tic plane crossings are shown in Fig. 3 as functions of the radial velocity V r . 
Case A gives rise to a wide range of ages between ~ 1 Myr and 210 Myr, and 
peculiar velocities between 110 km s -1 and 1500 km s -1 . Cases B and C, on 
the other hand, yield kinematic ages of at least ~ 55 Myr and ~ 125 Myr; and 
post-SN peculiar velocities of 90-270 km s _1 and 100-220 km s _1 , respectively. 

The post-SN orbital separation A and orbital eccentricity e at the times of 
Galactic disk crossings are obtained by numerical integration backwards in 
time of the equations governing the evolution of the orbit under the influence of 
gravitational radiation. The resulting post-SN orbital parameters range from 
A = 3.28 R Q and e = 0.274 when r kin = Myr to A = 3.36 R Q and e = 0.282 
when Thin = 210 Myr. These post-SN orbital parameters together with the 
post-SN peculiar velocities impose constraints on the pre-SN progenitor of 
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Fig. 4. Limits on the pre-SN progenitor of PSRB1534+12 and on the kick velocity 
imparted to the last-born NS. The left-hand panels show the constraints for case A, 
the middle panels for case B, and the right-hand panels for case C. Red regions 
correspond to solutions for which the pre-SN binary is detached, while blue regions 
indicate the additional solutions associated with mass-transferring systems. 

PSRB1534+12. These are derived in a similar way as for the progenitor of 
PSR J0737-3037 (except that the problem depends on only one free parameter, 
V r ). The results of our analysis are summarized in Fig. 4. 

Based on the observational constraints available in 2004, PSRB1534+12 may 
have been detached as well as semi-detached at the time the second NS was 
born. In order to illustrate this, the solutions for which no mass transfer takes 
place at the time of the second SN explosion are indicated by the red regions 
in Fig. 4, while the blue regions indicate the additional solutions that become 
accessible when the possibility of mass transfer is taken into account. For the 
latter solutions we adopt the same assumption as before that the system may 
survive the mass tra nsfer phase and form a DNS only if the mass ratio is 
smaller than 3.5 (e.g., Ilvanova et al.ll2003l ). 



For case A disk crossings, detached pre-SN binary configurations exist only for 
V r < — 250kms _1 . For larger radial velocities, the system is always undergoing 
mass transfer from the progenitor of the second-born NS to the first-born NS. 
In this analysis, the pre-SN mass of the helium star forming the second NS is 
constrained within 2.1 M & — 8 M Q . The lower limit again corresponds to the 
lowest mass for which a helium star is expected to form a NS instead of a 
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Fig. 5. Probability distribution functions of the magnitude of the kick velocity im- 
parted to PSRB1534+12's companion at the time of its formation, for cases A, B, 
and C. The left panels show the entire set of PDFs associated with all admissible 
Vr-values by means of a linear color scale that varies over yellow, green, orange, 
red, and blue with increasing PDF values, while the right panels show the PDFs 
associated with some specific Vr-values. For clarity, the curves in the right panels 
are offset from each other by an arbitrary amount. 

white dwarf, while the upper limit corresponds to the highest mass for which 
a helium star is expected to form a NS rather th an a black hole (see, e.g.. 
Fig. 1 in iBelczynski et al.l l2002bl and Table 16.4 in iTauris fc van den Heuvel 
20061 ). The divide between the red and blue regions at 4.7 M corresponds to 
the adopted critical mass ratio of 3.5 for dynamically stable mass transfer. 
The allowed mass range of pre-SN helium star masses is most constrained for 
\V r \ < 200 km s" 1 when 2.1 M < M <4M Q . Lower and upper limits on the 
pre-SN orbital separation are given by 2.4 R & and 4.3 _R . 



The behavior of the kick-velocity magnitude is similar as for PSR J0737-3039: 
the kick velocity increases with increasing absolute values of V r and, for a 
given radial velocity, it is generally constrained to an interval that is less than 
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Fig. 6. Probability distribution functions for the misalignment angle A between 
PSRB1534+12's spin axis and the post-SN orbital angular momentum axis, for 
cases A, B, and C (cf. Fig. 5 for more details). 

~ 600 km s -1 wide. Solutions for which the pre-SN helium star is overflowing 
its Roche lobe furthermore allow for somewhat higher kick velocities than 
solutions for which the helium star fits within its critical Roche lobe. These 
higher kick velocities are associated with the smaller orbital separations that 
become accessible for mass-transferring progenitors. When all possible pre-SN 
binary configurations are considered, the kick velocity magnitudes range from 
45kms~ 1 to 1350 km s -1 . When configurations for which the helium star is 
overflowing its Roche lobe are excluded, the range narrows to 230 km s -1 < 
Vk ^ 1350 km s -1 , so that the minimum required kick velocity is higher than 
when the possibility of Roche- lobe overflow is taken into account. The larger 
minimum kick velocity is required to compensate the larger effect of the mass 
lost from the system during the SN explosion: since the minimum Mo for 
these systems is 4.7 M Q , at least 57% of the total pre-SN mass is lost from the 
system. 

For case B disk crossings, almost all allowed pre-SN binary configurations 
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imply that the helium star progenitor of the second-born NS is overflowing 
its Roche lobe at the time of its SN explosion. The mass of the helium star 
is constrained to 2.1 M Q < M < 4.9 M , the pre-SN orbital separation to 
2.4 R Q < A < 4.3 R & , and the kick-velocity magnitude to 30kms _1 < Vk < 
475 km s -1 . If the helium star fits within its critical Roche lobe at the time of 
its SN explosion, the constraints become M ~ 4.8 M , A « 4.2 R Q , Vk ~ 
240 km s" 1 . 



Case C disk crossings, finally, yield no detached solutions for the progenitor 
of PSRB1534+12. The constraints for this case are similar to those found for 
case B, except that the mass Mo is always smaller than ~ 4 M and the kick 
velocity Vk is always smaller than 370 km s -1 . 



The constraints derived above may again be used to derive probability dis- 
tribution functions for the magnitude Vk of the kick velocity imparted to the 
second-born NS and for the misalignment angle A between the spin-axis of 
the first-born NS and the post-SN orbital angular momentum axis. We note 
that, unlike the PDFs derived for the double pulsar, the PDFs presented 
in this section do not incorporate probability distributions for V r obtained 
from theoretical population synthesis calculations, nor do they incorporate 
the probability of finding the system at its current place in the Galaxy. The 
resulting PDFs are presented in Figs. 5 and 6. For radial velocities of only 
a few lOOkms -1 , the kick- velocity distributions show two closely spaced and 
fairly evenly matched peaks between Vk — 50kms~ 1 and Vk — 250 km s -1 . For 
higher radial velocities, relevant only to case A, the peak(s) shift to larger kick 
velocities up to a maximum of ~ 1350 km s -1 . The tilt-angle distributions, on 
the other hand, favor misalignment angles below 30° when |V^| < 200kms _1 
and above below 100° when \V r \ > 600 km s -1 . For case A disk crossings, tilt 
angles close to A ~ 90° are furthermore strongly disfavored regardless of the 
value of the radial velocity. For case B and C disk crossings, tilt angles with 
non- vanishing probabilities are always smaller than 40°-50°. For kick velocities 
of less than 200kms _1 , the predictions in all three cases a r e in go od agreement 



with the measurement of A = 25.0 ± 3.8° by lStairs et al.l (120041 ) . 



We co nclude by noting that updated progenitor c onstraints bylThorsett et al. 
( 2005 ) incorporating the tilt angle measured by Stairs et al. ( 20041 ) exclude 



the possibility that the progenitor was detached right before the SN explosion 
forming the second NS. 



2.3 PSR B 1913+16: The Hulse-Taylor Binary Pulsar 



Besides the knowledge of the proper motion direction, PSRB1913+16 has the 
advantage that the misalignment angle between the pulsar's spin axis and the 
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pre-SN orbital angular momentum axis has been determined to be around 
~ 20°, correspon ding to prograde rotation, or ~ 160°, co r responding toret- 
rograde rotation ( Kramer|[l998l ; IWeisberg fc Tavlorl 12002 ). IWex et aD ( 2000l ) 
used this information to derive constraints on the mass of the second-born 
NS's direct progenitor, on the pre-SN orbital separation, and on the magni- 
tude and direction of the kick velocity imparted to the second-born NS at 
birth. In agreement with the then available numerical simulations of rapidly 
accreting neutron stars, the authors assumed that mass transfer from a helium 
star companion would cause the NS to collapse into a black hole and there- 
fore excluded Roche-lobe overflowing helium stars as vi able progen i tors o f 
the second-bo r n NS . How ever, more recent ca lculations by lDewi et al.l (120021 ). 
Ivanova et al.l (120031 ) , and iDewi fc Polsl (120031 ) show that NS binaries may sur- 
vive a helium-star mass-transfer phase, provided that the ratio of the helium 
star's mass to the neutron star's mass is not t oo ex treme (< 3.5). We therefore 
revise the constraints derived by IWex et al.l (120001 ) in the light of this new in- 
formation. The constraints on the tilt angle, for which we consider the values 
A = 18° ± 6° and A = 162° ± 6° are also imposed. 



We look for possible birth sites of PSRB1913+16 by tracing the system's 
motion in the Galaxy backwards in time up to a maximum age o f 80 MyiQD. 
as a fu nction of the unknown radial velocity V r . In agreement with lWex et al. 
(120001 ). we find that the system may have crossed the Galactic plane up to 



two times. The first Galactic plane crossing (case A) occurs at very young 
kinematic ages of ~ 2-4 Myr and gives rise to peculiar velocities in excess of 
~ 300kms _1 . The corresponding post-SN orbital parameters are A « 2.8 R Q 
and e ~ 0.618. The second Galactic plane crossing (case B), on the other hand, 
takes place at least ~ 55 Myr in the past and yields peculiar velocities of ~ 
230-440 km s _1 . The associated post-SN orbital separations and eccentricities 
range from A = 3.1 R Q and e = 0.646 to A = 3.2 R Q and e = 0.658. 



The constraints on the pre-SN parameter space accessible to the progenitor of 
PSRB1913+16 are shown in Fig. 7 as functions of the unknown radial velocity 
V r . As before, the red regions correspond to the solutions for which no mass 
transfer takes place at the time of the second SN explosion, while the blue 
regions indicate the additional solutions found when the possibility of mass 
transfer is taken into account. The constraints for detached pre-SN binary con- 



figura tions are in good agreement with the constraints derived by IWex et al 



(120001 ). It is clear however, that when the possibility of mass transfer is taken 
into account, the available parameter space becomes much less constrained. 



2 As for PSR J0737-3039 and PSRB1534+12, we use the spin-down age T b = 80 Myr 
instead of the characteristic age r c = 110 Myr as an upper limit for the age of the 
system. The maximum amount of orbital evolution that may hav e taken pla c e sinc e 



the formation of the DNS is therefore somewhat smaller than in Wex et al. I (|2000l h 
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Fig. 7. Limits on the pre-SN progenitor of PSRB1913+16 and on the kick veloc- 
ity imparted to the last-born NS. Red regions correspond to solutions for which 
the pre-SN binary is detached, while blue regions indicate the additional solutions 
associated with mass-transferring systems. 

For case A disk crossings, the range of radial velocities for which physically 
acceptable solutions exist is restricted to \V r \ < 500-600 km s -1 when A = 
18°, and to 100 kms" 1 < \V r \ < 1300 kms^ 1 when A = 162°. The mass of 
the second-born NS's direct progenitor is constrained to the interval between 
2.1 Mq and 8M , and the pre-SN orbital separation to the interval between 
1.1 Rq and 5.3 R Q for both considered values of the tilt angle A. The magnitude 
of the kick velocity varies from ~ 190 km s -1 to ~ 600 km s^ 1 when A = 18° 
and from ~ 580 km s -1 to ~ 2000 km s -1 when A = 162°. When the solutions 
are restricted to detached pre-SN binary configurations the limits become 
300 km s- 1 <V k < 600 km s" 1 and 700 km s" 1 < V k < 1250 km s" 1 for A = 18° 
and A = 162°, respectively. 

For case B disk crossings, the mass Mq of the second-born NS's direct progeni- 
tor and the pre-SN orbital separation Aq are constrained to the ranges of values 
given by 2.1 M < M < 8 M and 1.1 R Q < A < 5.3 R G when A = 18°, and 
to 2.1 M < M < 5M and 2.9 R Q < A < 5.3 R & when A = 162°. The 
magnitude of the kick velocity varies between ~ 190 km s -1 and ~ 530kms _1 
when A = 18°, and between ~ 550 km s -1 and ~ 850 km s -1 when A = 162°. 
When the solutions are restricted to those where no Roche-lobe overflow oc- 
curs at the time of the helium star's SN explosion, the minimum kick velocity 
associated with A = 18° increases slightly to approximately 280kms _1 , while 
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Fig. 8. Probability distribution functions of the magnitude of the kick velocity im- 
parted to PSRB1913+16's companion at the time of its formation, for both A = 18° 
and A = 162° (cf. Fig. 5 for more details). In the right-hand panels, the solid lines 
correspond to A = 18° and the dashed lines to A = 162°. 

the range of admissible kick velocities associated with A = 162° becomes very 
tightly constrained to ~ 640-680 km s _1 . 

The probability distribution functions for the magnitude of the kick velocity 
imparted to the second-born NS in PSRB1913+16 are shown in Fig 8 for both 
A = 18° and A = 162°. The distributions generally show a single peak at kick 
velocities which increase with increasing absolute values of V r . For A = 18°, 
the most probable kick velocity ranges from ~ 300 km s -1 to ~ 600kms _1 , 
while for A = 162° it ranges from ~ 600kms _1 to almost ~ 2000 km s -1 . 



3 Formation Channels for Double Compact Objects 



The inspiral and coalescence of double compact objects (DCO), such as NS- 
NS, BH-NS, and BH-BH binaries are some of the most promising candidate 
events for GW detection by current ground-based interferometers, like LIGO. 
The formation of double compact object populations has been studied by a 
number of dif ferent groups using pop ulation synthesis techniques. In one of 
these studies (IBelczynski et al.ll2002bl ) we focused on an investigation of the 
main formation channels, their origin, and their relative contributions. The de- 
tails of this study the synthes i s code used (StarTrack) and the conclusions are 
presented in IBelczynski et al. f |2002bl ). Here we merely summarize the findings 
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related to the various formation channels and their relative importance. Re- 
sults are described f or our standard model fr om that study (model A described 
in §2.1 and §2.2 of iBelczvnski et al.ll2002bh . 



We consider double compact objects with NS or BH (NS-NS, BH-NS or BH-BH 
binaries) with merger times shorter than 10 10 yr. In our standard evolutionary 
model, the population of coalescing DCO is dominated by NS-NS systems 
(61%), with a significant contribution by BH-BH binaries (30%), and a small 
contribution by BH-NS objects (9%). 

In Table 3, we present the most important formation channels of coalescing 
DCO, for our standard model. Formation channels of NS-NS, BH-NS and 
BH-BH binaries are marked by NSNS, BHNS and BHBH, respectively, they 
are listed in order of decreasing relative formation frequency (second column) 
with respect to the whole DCO coalescing population. The details of each 
evolutionary sequence, i.e., MT episodes and SN explosions are also given. 
Results were obtained based on the evolution of 3 x 10 7 primordial binaries. 



3.1 Populations of Double Neutron Stars 



Belczynski et al.l ( 12002bl ) identified a number of new NS-NS formation chan- 
nels. This is a result of tw o improvements in the implem entation of our popula- 
tion synthesis code, since IBelczvnski fc Kalogerall200ll . herea fter BK01. First, 
we ha ve replaced the approximate prescription suggested by lBethe fc Brown 
(119981 ) for the hyper-critical accreti on during CE phases, wi th a newly derived 
numerical solution (see Appendix of IBelczvnski et al.ll2002bl ). Second, we allow 
for hyper-critical CE evolution of low-mass helium stars with compact objects. 
In BK01, we allowed binaries with low-mass helium giants to evolve through 
DCE and standard SCE, but we had assumed that CE events of helium giants 
with compact ob jects lead to mergers, and possibly a gamma-ray burst (e.g., 
Fryer et al.lll998l ). However, due to the small mass of helium giant envelope 
at the onset of CE event (~ 1 — 1.5M ), we find that these systems survive 
the CE events, and form very tight NS-NS binaries. 



Double neutron stars are formed in various ways through more than 14 dif- 
ferent evolutionary channels identified in Table 3. The relative formation ef- 
ficiencies shown f or each channel are for the standard model A, described in 
Belczynski et al. fl2002bh . The entire population of coalescing NS-NS systems, 
may be divided into three main subgroups. 



Group I. This subpopulation consists of non-recycled NS-NS systems, first 
identified by BK01. These are systems in which none of the two NS ever had 
a chance of getting recycled through accretion. Our current results for the 
predicted formation rates and properties of the non-recycled NS-NS systems, 
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have not been affected by the two improvements discussed above. As shown 



Double Compact Object Formation Channels 

Formation Relative 
Channel Efficiency" 



Evolutionary History' 3 



NSNS:01 20.3 % NC:a^b, SN:a, HCE:b^a, HCE:b^a, SN:b 

NSNS:02 10.8 % NC:a^b, SCE:b^a, NC:a^b, SN:a, HCE:b^a, SN:b 

NSNS:03 5.5 % SCE:a^b, SN:a, HCE:b^a, HCE:b^a, SN:b 

NSNS:04 4.0 % NC:a^b, SCE:b^a, SCE:b^a, SN:b, HCE:a^b, SN:a 

NSNS:05 3.2 % DCE:a^b, SCE:a^b, SN:a, HCE:b^a, SN:b 

NSNS:06 2.5 % SCE:a^b, SCE:b^a, NC:a^b, SN:a, HCE:b^a, SN:b 

NSNS:07 2.2 % NC:a-^b, NC:a-*b, SN:a, HCE:b^a, HCE:b^a, SN:b 

NSNS:08 2.0 % NC:a^b, DCE:b^a, SN:a, HCE:b^a, SN:b 

NSNS:09 2.0 % DCE:a^b, DCE:a^b, SN:a, SN:b 

NSNS:10 1.6 % NC:a^b, SCE:b^a, SN:b, HCE:a^b, SN:a 

NSNS:11 1.5 % NC:a^b, SCE:b^a, DCE:b^a, SN:a, SN:b 

NSNS:12 1.5 % NC:a^b, SCE:b^a, DCE:a^b, SN:a, SN:b 

NSNS:13 1.0 % DCE:a^b, SN:a, HCE:b^a, SN:b 

NSNS:14 3.0 % all other 



BHNS:01 
BHNS:02 
BHNS:03 
BHNS:04 



4.5 % 

1.6 % 
1.3 % 
2.0 % 



NC:a^b, SN:a, HCE:b^a, SN:b 
NC:a^b, SCE:b^a, SN:a, SN:b 
SCE:a^b, SN:a, HCE:b^a, NC:b^a, SN:b 
all other 



BHBH:01 
BHBH:02 
BHBH:03 



17.7 % 
10.5 % 
1.4% 



NC:a^b, SN:a, HCE:b->a, SN:b 
NC:a-^b, SCE:b^a, SN:a, SN:b 
all other 



"Normalized to the total DCO population. 

^Sequences of different evolutionary phases for the primary (a) and the secondary 
(b): non-conservative MT (NC), single common envelope (SCE), double common 
envelope (DCE), common envelope with hyper-critical accretion (HCE), supernova 
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explosion/core-collapse event (SN). Arrows mark direction of MT episodes. 



in Table 3, the non-recycled NS-NS systems are formed via the NSNS:09, 
NSNS:11 and NSNS:12 channels, which involve DCE of two low-mass helium 
giants that were already allowed in the earlier version of StarTrack. 

The unique qualitative characteristic of this NS-NS formation path is that both 
NS have avoided recycling. The NS progenitors have lost both their hydrogen 
and helium envelopes prior to the two supernovae, so no accretion from winds 
or Roche-lobe overflow is possible after NS formation. Consequently, these 
systems are detectable as radio pulsars only for a time (~ 10 6 yr) much shorter 
than recycled NS-NS pulsar lifetimes (~ 10 s — 10 10 yr in the observed sample). 
Such short lifetimes are of course consistent with the number of NS-NS binaries 
detected so far and the absence of any non-recycled pulsars among them. 

We note that the identification of the formation path for non-recycled NS- 
NS binaries stems entirely from accounting for the evolution of helium stars 
and for the possibility of double CE phases, both of whi ch have typically 



been ignored in previous calculations (with the exception of iFryer et al.lll998l . 
where, however, such events were assumed to lead to mergers). 

Group II. This subpopulation consists of tight, short lived binaries with one 
recycled pulsar. Their merger times are typically ~ 1 Myr or even smaller (see 
§3.4.5). As shown in Table 3, these new dominant NS-NS systems are formed 
via the NSNS:01-08, NSNS:10 and NSNS:13 channels, with the common char- 
acteristic that the last binary interaction is a hyper-critical CE of a low-mass 
helium giant and the first-born NS. 



In iBelczynski et al.l ( j2002al ) we describe in detail the formation of a typical 



NS-NS binary of group II. The channel identified as the most effic ient for NS 



NS for mation (NSNS:01) corresponds to the "standard" channel of IFryer et al. 



( 119981 ). The only difference is an extra CE event which originates from allowing 
for helium star evolution and without a priori assumptions about the CE 
outcome. The second most dominant channel, involving two consecutive MT 
episodes and then two SN explosions, closely resembles our channels: NSNS:02, 
NSNS:04, NSNS:06, NSNS:10, NSNS:11, NSNS:12. The only difference again 
remains an extra MT episode from evolved, Roche-lobe-filling helium stars. 

The most dramatic effect of the binary evolution updates is reflected in the ex- 
istence of a whole new population of coalescing NS-NS stars formed in Group 
II. In our standard model these channels contribute 50% of the DCO popula- 
tion, and their common characteristic is that the last binary interaction is a 
hyper-critical CE of a low-mass helium giant and the first-born NST 3 !. It turns 



3 A more careful treatment of the response of helium stars to mass loss leads to a 
small reduction of this percentage, but the population forming through this channel 
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out that the majority of these systems survive the HCE event and form tight 
NS-NS binaries. Had we not taken into account the radial expansion of low- 
mass helium-rich giants, the progenitors of this dominant NS-NS population 
would have evolved without any further MT. Most of them would have still 
formed NS-NS systems, although not as tight as after this last CE episode. 
We have actually examined this alternative and found that about half of them 
would have formed binaries with merger lifetimes longer than 10 10 yr. 

Group III. This subpopulation consists of all the other NS-NS systems formed , 
through more or less classical channels (IBhattacharya fc van den Heuvellll99ll ). 
The formation path denoted NSNS:14 corresponds to what is usually consid- 
ered to be the "standard" NS-NS formation channel (IBhattacharya fc van den Heuvel 

199ll ). Since we account for hyper-critical accretion in CE, the formation rate is 

decre ased be cause some NS (bu t not all, as assumed by lPortegies Zwart fc Yungelson 
19981 and by lFryer et al.lll998l ) collapse to BH. Furthermore our treatment of 
the hyper-critical accretion typically leads to tighter post-CE systems, caus- 
ing more binaries to merge in CE events, and thus decreases the number of 
possible NS-NS progenitors. 



Brown! (119951 ) proposed a NS-NS formation channel where the progenitor con- 
sists of two massive stars with nearly equal masses (within 4%) and the first 
mass transfer episode occurs when both stars are on the giant branch and 
have convective envelopes; if the binary survives the ensuing common enve- 
lope phase, both envelopes are ejected (double common envelope; DCE) and 
two helium stars are exposed which explode as supernovae and a NS-NS is 
formed. Brown was motivated by the almost equal mass measurements in NS- 
NS binaries and by adopting a low maximum NS mass (1.5 M ; associated 
with soft equations of state). In the proposed channel the progenitors are 
of almost equal mass and none of the NS ever experiences a common enve- 
lope phase, avoids collap sing into a BH. Some of the channels identified by 
Belczynski etall ( l2002bf ) (Table 3; specifically NSNS:09,11,12) are similar to 
the iBrownl (119951 ) channel. However for the standard model they contribute 
only 8% of the total NS-NS population because a more widely accepted max- 
imum neutron star mass of 3 M is adopted and NS do not collapse into BH 
even if they experience hypercritical accretion in common envelopes (HCE). 
Recently Pinsonneault & Stanek (2006) reported on a population of "twin" 
binaries with component masses within 5%; this population about 50% of the 
observed sample, but of course it is also favored by selection ef fects and the true 



contrib ution of "twins" could be closer to 25% instead of 50%. lBelczynski et al 



( I2002bl ) did also consider a model with a maximum NS mass of 1.5 M (model 
D2); as expected the NS-NS rate is decreased significantly especially in com- 
parison to BH-NS, since most of the standard NS-NS channels in Table 3 
actually form BH-NS (the rate ratio shifts from 6.5 for the standard model A 



is still significant, as discussed in llvanova et all (|2003h and iBelczvnski et alj (2006 
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to 0.25 for model D2; see Table 4 in iBelczynski et al.ll2002bl ). We note that a 
more recent population study of the Brown (1995) channel presented by Dewi 
et al. (2006) leads to NS-NS rates consistent with that from model D2. 



3.2 Populations of Black Hole Binaries 



In general, BH-NS and BH-BH binaries are formed through just a few dis- 
tinct channels, with a moderate number of MT events (2-3), in contrast to 
our findings for NS-NS systems. Helium star evolution, radial expansion and 
CE phases are much less important for their formation. The reason is that 
for most of these progenitors the first-born compact object is massive enough 
that they do not expand to large radii nor they lead to possible CE evolution 
(see channel BHNS:03 in Table 3). Instead these DCO form most efficiently 
through channels t hat closely resemble those NS-NS con v entionally though t 
to be "standard" (IBhattacharya fc van den Heuvell Il99ll ; iFryer et al.l Il998l ): 
evolution is initiated with a phase of non-conservative mass transfer and fol- 
lowed either by a CE phase or the formation of the first compact object (see 
BHNS:01, BHNS:02, and BHBH:01, BHBH:02). 



4 Merger Rates of Double Neutron Stars: Observed Pulsar Sample 



4-1 Introduction 



Double Neutron Stars that will merge within a Hubble time are one of the 
prime targets for ground-based gravitational- wave (GW) interferometers such 
as GEO600, TAMA, VIRGO, and LIGO. Event rates of the DNS inspiral 
searches by these detectors can be inferred using the rate estimates with an 
extrapolation out to the maximum detection distances for any detector under 
consideration. Before 2003, the Galactic DNS merger rate had been estimated 



yr" 



sec 



between ~ 10 7 — 10" 

At that time, th ere were only two sy stems available for empirical stu dies, 



Kalogera et al.ll200ll and references therein). 



PSRs B1913+16 flHulse fc Tavlorlll975t ) and B1534+12 flWolszczanlll99lh . We 
calculated th e PDF of Galact ic DNS merger rates, P(R), based on these two 



systems (see 



Kim et al. 2003 . hereafter KKL, for further discussiorF"!) . Soon 
after the discovery of the highly relativistic system PSR J0737— 3039 , we were 
able to revise P(R) including PSR J0737— 3039 in collaboration with the ob- 



4 We note that a web interface that allows the calculation of P(R) 
for a wide range of pulsar surveys and pulsars is publicly available at 



ht tp : / / www .astro . northwestern . edu/ ~ciel/ gppg_main . html 
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servation teanf^lf Burgay et al. 2003 ; Kalogera et al. 2004a ). The discovery of 
J0737— 3039, resulted in a significant increase in the estimated Galactic DNS 
merger rate by a factor of ~6. This implies a boost in event rates for DNS 
searches for GW interferometers. Here, we present our recent results on: (i) 
the PDF of Galactic DNS merger rate estimates with updated observations; 

(ii) the global PDF of rate estimates considering the systematic uncertainties; 

(iii) constraints on upper limits for rate estimates based on the observed su- 
pernova rate incorporated with our theoretical understanding of the SN-DNS 
relation; (iv) the approximate contribution of J1756— 2251 to the Galactic 
DNS merger rates and uncertainties in our assumptions on the efficiency of 
the acceleration searches for Parkes multibeam pulsar survey (PMPS). Fi- 
nally, we dis cuss implications of the most recently discovered pulsar binary 
J1906+0746 dLorimer et alJl2006t ). 



4-2 The Galactic DNS merger rate 



Here, we describe the main components of the calculation of the combined 
P(R) considering the three observed DNS systems in the Galactic disk. Th e 
details of these calculations are discussed by KKL and (IKalogera et al.ll2004al ). 
The merger rate of a given population % can be defined by 



Ra 



Npsr 

Tlife 



fl 



(3) 



where iVpsR,; represents the number of pulsars in our Galaxy with pulse and 
orbital characteristics similar to an observed sample i (e.g. PSR J0739— 3039) 
and fhi is a correction factor to take into account pulsar beaming (typically 
~ 6 j_j; Tiif ei i is the lifetime of system i based on its observed properties. In or- 
der to calculate iVp SR i , we distribute a large population of pulsars all similar 
to the system i, in a model galaxy assuming spatial and luminosity distribu- 
tions. Since pulsar luminosities are drawn from a distribution, the observed 
flux and estimated distance for the DNS system are not relevant in our cal- 
culation. Moreover, the selection effects for faint pulsars are implicitly taken 
into account. Once we generate a pulsar population with a size iVp SR , we can 
then calculate the number of pulsars detected (iV det ) by large-scale pulsar sur- 
veys. We repeat the survey simulations with a detailed modeling of selection 
effects for observed DNS systems. For a fixed iVpsR, we find that A^et follows 
a Poisson distribution, P(A^ e t; N^et), where iVdet is a mean value of iVdet for a 
given population size (Npsr). We require iVdet = 1, i.e., we consider only one 

5 Only the millisecond component (J0737-3039A) is considered in our calculation. 
For instance, the current age of the system is derived from pulsar A. 

6 This is based on polarimetry measurements of B1913+16 and B1534+12. 
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observed system, and calculate the best-fit value of iVdet- With a wide range 
of iVp SR , we find N det = cN PSR as expected where c is a constant. Applying 
Bayes' theorem to these results, we calculate a P(iVpgpJ, a PDF for the pop- 
ulation size of a given system i knowing that there is one observation. The 
calculation P(R) from P(Npsr) is straightforwad, and the full derivatio n and 
the analytical formula can be found in Appendix A of iKim et al.l ( 120041 ) . 



The lifetime of a DNS system is defin ed by 7iif P = t s h + TWg, where r sd is 



the spin-down age of a recycled pulsar (lArzoumanian et al.lll999[) and T mr a is 



the r emaining lifetime until the two neutron stars merge (IPeters & Mathews 
19631 ) . Based on the most re cent observations , we estimate the lifetime of 
J0737-3039 to be ~230Myr flLvne et al.ll2004h . This is the shortest among 
the three observed systems. The beaming correction factor / b is defined as 
the inverse of the fractional solid angle subtended by the pulsar beam. Its 
calculation requires d etailed geometrical information on the beam. Following 
Kalogera et al.l (120011 ) and updating the values with recent ob servations, we 
calcualte f h ~5.7 for PSR B1913+ 16 flWeisberg fc Tavlorl l2002h and ~6.0 for 
PSR B1534+12 flStairs et al.ll2004 ) Without good knowledge of the geometry 
of J0737— 3039, we assume /b,jo737 to be the average value of the other two 
systems (~ 5.9). 



In Fig. 9, we show P{R) for a reference model (Model 6 in KKL). We obtain the 
most likely value of R ~ 71 Myr -1 , larger by a factor of ~ 5.5 than the rate 
estimated before the discovery of J0737— 3039. The increase factor is found 
similar for all pulsar population models we examined. The increased merger 
rates imply a boost in the inferred detection rate of GWs from DNS inspirals 
for ground-based GW interferometers such as LIGO. In order to calculate 
the detection rate (D), we assume a homogeneous distribution of galaxies 
in nearby Universe and a spherical symmetry in detector sensitivity. Then, 
we can write D = i? gal x iV gal , where N gal is the number of galaxies in the 
detection volumne (Vd e t)- We calculate the number density of galaxie s derived 



by the obse rved blue luminosity de nsity, (n ga i = 1.25x 10 _2 Myr _1 (see lPhinney 
( 1991 ) and Kalogera et al.l ( 2001 ) for more details). The detection volume of 
LIGO can then be defined as a sphere for a given detection distance (20 Mpc 
and 350 Mpc for the inital and advanced LIGO, respectively), and the number 
of galaxies within Vdet is simply n ga i x V^et- F° r our reference model, we find 
that the most probable event rates are about 1 per 30 yrs and 1 per 2 days, for 
initial and advanced LIGO, respectively. At the 95% confidence interval, the 
most optimistic predictions for the reference model are 1 event per 9 yrs and 
1.6 events per day for initial and adva nced LIGO, respectively. More details 
can be found in iKalogera et al.l (j2004al ). 



As shown in Fig. 9, the Galactic DNS merger rate is dominated by PSR 
J0737-3 039. We note that th e current age of 30-70Myr for J0737-3039, sug- 
gested by lLorimer et al.l (120051 ). implies a even shorter lifetime (ri ife ~ 115—155 
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Fig. 9. P(R) is shown on a log scale. The thick solid line is the Galactic rate estimate 
overlapped with results for individual observed systems (dashed lines). Dotted lines 
indicate confidence intervals for the rate estimates. The same results are shown on 
a linear scale in the small inset. All results are from our reference model. 

Myr) knowing that the estimated merger timescale of this system is ~85 Myr. 
Based on their results, we find the most likely value of R for the reference 
model is ~ 90 - 110 Myr" 1 . 

The beaming correction for J0737— 3039 is not yet constrained and we assume 
MSPs discovered in DNS systems are not very different. As a conservative 
lower limit, without any beaming corrections for all observed systems for a 
reference model, we obtain ~ Myr" 1 with a 95% confidence interval. The 
corresponding detection rates for initial and advanced LIGO are 5lj 2 x 10" 3 
yr" 1 and 27iJi y r_1 > respectively. Only when the axes geometry of J0737- 
3039 becomes available, we will be able to constrain the uncertainties of the 
beaming fraction, and in turn, the rate estimates. 
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4-3 Global probability distribution of the rate estimates 



In KKL, we showed that empirical DNS merger rates are strongly dependent 
on the assumed luminosity distribution function for pulsars, but not on the 
pulsar spatial distribution. Therefore, we can consider only the rate depen- 
dence on the pulsar luminosity function for simplicity. Here, we describe how 
we can incorporate the systematic uncertainties from these models and cal- 
culate, P g (R), a global PDF of rate estimates. Note that the results available 
on prior functions for pulsar luminosity distribution are currently out of date. 
Specific quantitative results could change when constraints on the luminosity 
function are derived from the current pulsar sample. 

We assume a power-law luminosity distribution for a radio pulsar luminosity 
function f(L). This function is defined by two parameters: the cut-off lumi- 
nosity L m in and power-index p. We assume prior distributions for these two 
parameters an d calculate P S (R). We fit t he marginalised likelihood of L m i n and 



p presented by lCordes &: Chernofil (119971 ) and obtain the following analytic for- 
mulae for prior functions, i.e. f(L min ) and g{p): /(-^min) = tto+^i^mm+c^min 
and g(p) = lO /3o+/3lP+/32P , where cii and $ (i = 0, 1,2) are coefficients we ob- 
tain from the least-square fits and the functions are defined over the inter- 
vals Lmin = [0.0, 1 . 71 mJ y kpc 2 and p = [1.4, 2.6]. We note that, although 
Cordes fc Chernofil (119971 ) obtained f(L min ) over L min ~ [0.3, 2] mJy kpc 2 



centered at 1.1 mJy kpc 2 , we consider f(L min ) with a peak at ~ 0.8 mJy 
kpc 2 consideri ng the discoveries of faint pulsars with L 140 o below 1 mJy kpc 2 



(jCamiloll2003l ). Next we calculate P g {R) as follows: 



P S (R) = J dp J dL min P(R)f(L min )g(p). (4) 



In Fig. 10, we show the distributions of L min and p adopted (top panels) 
and the resulting global distribution of Galactic DNS merger rate estimates 
(bottom panel). We find the peak value of P g (R) at only around 13Myr _1 . 
We note that this is a factor ~ 5.5 smaller than the result from our reference 
model (R ~71Myr _1 ). At the 95% confidence interval, we find that the global 
Galactic DNS merger rates lie in the range ~ 1-145 Myr -1 . These imply LIGO 
event rates in the range ~ (0.4 — 60) x 10~ 3 yr -1 (initial) and ~ 2 — 330 yr -1 
(advanced). Since 1997, the number of known millisecond pulsars has more 
than doubled, and therefore, constraints on L min and p and their PDFs based 
on the most up-to-date pulsar sample are urgently needed. 
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Fig. 10. The global Pg(R) on a linear scale (lower panel) and the assumed intrinsic 
distributions for L m \ n and p (upper panels). Dotted lines represent the lower (SNl) 
and upper (SN\j) bounds on the observed SN Ib/c rate scaled by 5% of the observed 
SN Ib/c rates, 600 - 1600 Myr" 1 (see text). 



4-4 Upper Limit of DNS Merger Rate Estimates. Constraints from Type Ib/c 
Supernovae Rates 



According to the standard binary evolution scenario, the progenitor of the sec- 
ond neutron star is expected to form during a Type Ib/c supernova. Therefore 
the empirical estimates for the Type Ib/c SN rate in our Gala xy can be used to 



set up per limits on the DNS merger rate estimates. Based on lCappellaro et al. 
(119991 ) . we adopt i?sNib/c — 1100 ± 500 Myr -1 considering Sbc-Sd galaxies, a 
Hubble constant H = 71km/s/Mpc and the blue luminosity of our Galaxy 
Lb, gai = 9 x 10 9 L sun . 
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We calculate the fraction of SN Ib/c actual ly involved in the formation o f 
DNS with a binary evolution code StarTrack feelczvnski et al.ll2002bl . l2006bh 



and estimate the rate ratio: 7 = (i?DNs/-RsNib/c) x 100 < 5%. Motivated by 
this result, we adopt the empirical i?sNib/c assuming 7 ~ 5% and compare the 
value with the global PDF (Fig. 10). 

We note that our most optimistic DNS merger rate is R = 189^ig6 Myr -1 at a 
95% confidence interval (Model 15 in KKL). We obtain 7 ~ 80% with respect 
to the center value of the empirical SN Type Ib/c rate (1100 Myr -1 ) and the 
upper limit of R at the 95% confidence interval (189 + 691 = 880 Myr -1 ). 
This corresponds to 7 ~ 13% for a SN Type II rate, which is factor 6.1 larger 
than that of SN Type Ib/c. In both cases, the most optimistic model is lower 
than the current empirical supernova rate estimates, but not really consistent 
with the results of population synthesis calculations. If we consider an upper 
limit of R at the 95% confidence interval from the global PDF, we obtain 
7 ~ 13% and 2% for the center value of SN Type Ib/c (1100 Myr" 1 ) and II, 
respectively. 



4-5 Implications of New Discoveries to the Galactic DNS merger rate esti- 
mates 



Recently. iFaulkner et al.l (120051 ) discovered PSR J1756— 2251, the 4th merging 
DNS in the Galactic disk from the Parkes Multibeam Pulsar Survey (PMPS). 
The standard Fourier method failed to find this pulsar and they reanalysed the 
PMPS data with an acceleration search (or a 'stack search' as described in their 
paper). In order to calculate the merger rate including PSR J1756— 2251, a 
detailed simulation is necessary to calculate the effect of the acceleration search 
with the PMPS. However, the approximate contribution of PSR J1756— 2251 
to the Galactic merger rate can be easily obtained. We find the total rate 
increases by only ~ 4% due to the new discovery. This is expected because PSR 
J1756— 2251 can be identified as a member of the B1913+16-like population, 
which has already been taken into account in the calculation. Only future 
detections of pulsars from a significantly different population (compared to 
the known systems), or from the most relativistic systems, will result in a 
non-trivial contribution to the rate estimates. 

Finally, we note the implications of J1906+0746 on the pulsar binary merger 
rates. This system has drawn attention du e to the extremely young age of 



the pulsar (characteristic age of ~ 112 kyr; lLorimer et al.ll2006l ). If the com- 
panion is another neutron star, J1906+0746 would be the first discovery of 
a non-recycled component in a DNS system. Currently, the nature of the 
companion is totally unknown, and it can be either a light neutron star, or 
heavy (O-Mg-Ne) white dwarf. Assuming J1906+0746 is a DNS system, we 
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calculate its contribution to the Galactic merger rate. Because of its short 
lifetime (~82Myr), J1906+0746 can increase the Galactic DNS merger rate 
by about a factor 2. This implies that the current estimated DNS merger 
rate including J0737— 3039 can still be d oubled! If J1906+0746 is an ecc en- 
tric NS-W D system, such as Jl 141-6545 (Ivan Kerkwijk fc Kulkarnilll999l ) or 
B2303+46 flStokes et alJll985h . it will be as important as J1141-6545, which 
is currently dominate the birthrate of eccentric NS-WD binaries. 



5 Merger Rates of Double Compact Objects: Population Synthesis 



Lacking a sample of DCOs containing a black hole, the only route to BH-BH 
merger rates is via population synthesis models. These involve a Monte Carlo 
exploration of the likely life histories of binary stars, given statistics governing 
the initial conditions for binaries and a method for following the behavior of 



single and binary stars (see, e.g., iBelczvnski et al.ll2002bl ; iHurley et al.l 12002 



Portegies Zwart Verbuntl Il996l ; iFryer et al.l Il998l . and references therein) . 



Unfortunately, our understanding of the evolution of single and binary stars 
is incomplete, and we parameterize the associated uncertainties with a great 
many parameters (~ 30), many of which can cause the predicted DCO merger 
rates to vary by more than one or two orders of magnitude when varied inde- 
pendently through their plausible range. To arrive at more definitive answers 
for DCO merger rates, we must substantially reduce our uncertainty in the pa- 
rameters that enter into population synthesis calculations through comparison 
with observations when possible. 



The simplest and most direct way to constrain the parameters of a given pop- 
ulation synthesis code is t o compare several of i ts man y predictions against 
observations. For example, lO'Shaughnessy et al.l (l2005bl ) required their popu- 
lation synthesis models be consistent with the empirically estimated formation 
rates derived from the four known Galactic NS— NS binaries which are tight 
enough to merge through the emission of gravitational waves within 10 Gyr. 
Here, we require our population synthesis models to be consistent (modulo 
experimental error) with six observationally determined rates: (i) the forma- 
tion rate implied by the known Galactic merging NS— NS binaries, mentioned 
above; (ii) the formation rate implied by the known Galactic NS— NS bina- 
ries which do not merge within 10 Gyr (henceforth denoted "wide" NS— NS 
binaries); (iii,iv) the formation rate implied by the sample of merging and ec- 
centric WD-NS binaries; and finally (v,vi) the type II and type Ib/c SNe rates. 
Further, we use the set of models consistent with these constraints to revise 
our population-synthesis-based expectations for various DCO merger rates, 
assuming no prior information, so all population synthesis model parameters 
consistent with our constraints are treated equally. 
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5.1 Observations of DCOs 



Seven NS— NS binaries and four WD-NS binaries (with relatively massive 
WDs) have been discovered so far in the galactic disk, using pulsar surveys 
with well- understood selection effects. [We are very s pecifically not includ - 
ing the recently-discovered binary PSR J1906+0746 (ILorimer et al.l 120061 ). 
because the companion cannot be definitively classified as a WD or NS at 
present. We also omit PSR J2127+11C found in the globular cluster M15, 
since its formation is thought to be dynamical and not due to isolated binary 
evolution in the Galactic field. 



5.2 Methodology I: Preferred population model 



As already discussed, in the context of NS— NS KKL developed a statistical 
method to estimate the likelihood of DCO formation rates, given an observed 
sample of DCOs in which one member is a pulsar, designed to account for 
the small number of known systems and their associated uncertainties. In this 
section, we present results only for our reference pulsar luminosity distribution 
model, corresponding a power law luminosity distribution with negative slope, 
index p = 2, and minimum luminosity L m i n = 0.3 mJy kpc 2 [model 6 of KKL; 
as discussed therein, this model better accounts for more recent observations of 
faint pulsars] ; in the following section, we describe our our predictions change 
when systematic uncertainties in p and L min are incorporated into V . 

Finally, for each class of binary pulsars we define symmetric 95% confidence 
intervals: the upper and lower rate limits 1Z W} ± satisfy 

7Zw,— oo 

J dKV w (K) = J dUV w {TZ) = 0.025 . (5) 

7^,4- 



This confidence-interval convention is different from with the customary choice 
presented in KKL. 

In all plots that follow, rate probability distributions are represented using a 
logarithmic scale for 71; thus, instead of plotting V, all plots instead show 



p(logft) = V(K)1l\nl0 . 
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5.2.1 WD-NS binaries 



Four WD-NS binaries with relatively massive WDs have been discovered in 
the galactic disk: PSRs J0751 + 1807, J1757-5322, J1141-6545, and B2303+46. 
While all four binaries may be applied to a net WD-NS formation rate esti- 
mate, the sample manifestly contains relics of distinctly different evolutionary 
channels; for example, while J0751+1807 and J1757-5332 have evidently been 
strongly circularized and spun up by a recent mass transfer e pisode, J1141- 



6545 cannot have been (IKaspi et al.l I200CH : iBailes et al.l 120031 ). Ideally, pop- 



ulation synthesis must produce distributions of WD-NS binaries consistent 
with both the observed orbital parameters and spins; however, the present 
sparse sample does not allow a reliable nonparametric estimate of the distri- 
bution of WD-NS binary parameters. Instead, as first step towards applying 
constraints based on binary parameter distributions, we subdivide these five 
binaries into two overlapping classes: merging binaries, denoted WD-NS (m), 
which will merge through the emission of gravitational waves within 10 Gyr; 
and eccentric binaries, denoted WD-NS (e), which have significant (e > 0.1) 
eccentricity at present. The rate estimate derived for b oth classes is dominated 



by J1141-6545 flKim et al.l 12004 : iKalogera et al.ll2005h . 



5.2.2 NS-NS binaries 

Seven NS— NS binaries have been discovered so far in the Galactic disk. Four 
of the known systems will have merged within 10 Gyr (i.e., "merging" bi- 
naries: PSRs J0737-3039, B1913+16, B1534+12, and J1756-2251) and three 
are wide with much longer merger times (PSRs J1811-1736, J1518+4904, 
and J182 9+2456). PSR J1756- 2251 was discovered recently with acceleration 



searches (IFaulkner et al.ll2005l ). As already discussed, the rate increase is es- 
timated to be smaller than 4%. For this reason, we omit it when estimating 
NS-NS merger rates. 

The observed NS-NS population naturally subdivides into two distinct classes, 
depending on whether they merge due to the emission of gravitational waves 
within 10 Gyr: merging NS— NS binaries, denoted NSNS(m), and wide NS— NS 
binaries, denoted NSNS(vw). 



5.2.3 Recycling, selection effects, and the lack of wide NS-NS binaries 

We find that the confidence intervals for wide and merging NS-NS binaries are 
almost an order of magnitude from overlapping. On the contrary, population 
synthesis simulations produce merging and wide binaries at a roughly equal 
(and always highly correlated) rates. Since our rate estimation technique auto- 
matically compensates for the most obvious selection effects (e.g., orientation 
and detectable lifetime), two unbiased samples of wide and merging NS-NS 
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binaries should arrive at nearly the same prediction for the NS-NS formation 
rate. 



O'Shaughnessy et al.l (j2005bl ) explained this significant discrepancy by a se- 



lection effect: evolutionary tracks leading to wide NS-NS binaries should be 
less likely to recycle the first-born pulsar. The observed sample of wide NS-NS 
binaries [NSNS (vw)] represents a much sma ller subset of recycled pulsars. As 
summarized in lO'Shaughnessy et al.l (j2006al ). we assume any wide NS-NS bi- 
nary which in its past underwent any conventional mass transfer will recycle 
one of its neutron stars to a pulsar. This condition is likely to over-estimate 
the true NSNS(vw) formation rate. We note that neutron-star recycling can 
also possibly happen during dynamically unstable mass-transfer phases (com- 
mon envelopes), but the vast majority of wide NS-NS do not evolve through 
such a phase. 



5.2.4 Pulsar Population Model Uncertainties 



As noted in KKL and subsequent papers, our reconstruction of the pulsar pop- 
ulation (i.e., Npsr) relies upon our understanding of pulsar survey selection 
effects and thus on the underlying pulsar luminosity dis tribution. This dis- 
tribu tion can be well-constrained experimentally (see, e.g. ICordes fc Chernofi 
19971 ) . though these constraints do not yet incorpo rate recent faint pulsar 



discoveries such as J1124-5916 flCamilo et al.l l2002h . Nonetheless, different 



observationally-consistent distributions imply significantly differen t distribu- 
tions , with maximum-likelihood rates differing by factors of order 10 (IKim et al 
20031 ) . Since the constraint intervals discussed above assume the preferred pul- 
sar luminosity distribution model - a power-law pulsar luminosity distribution 
with negative slope p = 2 and minimum cutoff luminosity L min = 0.3 mJy kpc 2 
- they do not incorporate any uncertainty in the pulsar luminosity function. 



The infrastructure needed to incorporate uncertainties in the pulsar luminosity 
function has been presented in §4.2. However, out-of-date constraints on the 
pulsar luminosity function allow implausibly high minimum pulsar luminosi- 
ties L min . A high minimum pulsar luminosity implies fewer merging pulsars 
have been missed by surveys. Thus, these out-of-date constraints on pulsar lu- 
minosity functions permit models consistent with substantially lower merger 
rates than now seem likely, given the dis covery of fa i nt pul sars. In other words, 
if we use the infrastructure presented in IKim et al.l (120061 ) to marginalize over 
L min and p generate a net distribution function for the merger rate, then the 
95% confidence intervals associated with that net distribution would have a 
spuriously small lower bound, entirely because the pulsar population model 
permits large L m i n . Therefore, we present results based only on our preferred 
luminosity fu nction and do not include the out-of-date luminosity function 
constraints by ICordes fc Chernofil (119971 ) and the related net rate distribution 
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provided by iKim et al.l (120061 ). 



5.3 Observations of supernovae 



Type Ib/c and II supernovae occur extremely rarely near the Milky Way. 
While historical data contains several observations of and even surveys for su- 
pernovae, the selection effects in these long-duration heterogeneous data sets 
have made their interpretation difficult (Cappellaro 2005, private communica- 
tion Vjjijjfis_rjaj^^ uncertainties via Table 
4 of bappellaro et all ril999h . ICappellaro et all <jl999h present their results m 
terms of a number of supernovae per century per 10 10 blue solar luminosities 
(L 0v b); to convert their results to rates per Milky Way equivalent galaxy, we 
assume a Milky Way blue luminosity of 2 x 10 10 L QjB = 0.9 x W W L^ relativ e 
to the solar blue (L^r ) and t otal ( L@) luminosities (see, e.g., iKalogera et e 
( 2001 ). Phinney ( 199ll ). Cox ( 2000l ) and references therein). Using their esti- 



mates for the most-likely supernovae rates and for the la errors in those rates 
for Milky Way-like galaxies (i.e., SOa-Sb), we arrive at the 2a logarithmic 
confidence intervals used in our constraints. 



Though fairly accurate studies exist of the high-redshift supernova rate (e.g., 
Cappellaro et al.ll2005l ). they have little relevance to the present-day Milky 
Way. Several surveys have also attempted to determine the supernova rate 
in the Milky Way by a variety of indirect methods, such as statistics of su- 
per nova remnants (highly unreliable d ue to challenging selection effects; see, 



e.g.. Ivan den Bergh fc Tammannlll99ll ) and direct observation of r adioisotope- 



prod uced backgrounds (e.g., decay from b Al, as described in iDiehl et al 
20061 ). Taken independ ently, these methods ha ve greater uncertainties than 
the historical studies of ICappellaro et al.l (119991 ) . 



5.4 Population synthesis predictions 



5.4-1 Star Track population synthesis code 



We estimate formation and merger rates for several classe s of double com- 
pact ob jects using the StarTrack code first developed by iBelczynski et al 
(I2002bf) and recently signifi cantly updated and tested as described in detail 
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Belczynski et al.l ( I2006bl ). This updated code predicts so mewhat different 



double compact object properties than the version used in IBelczynski et al. 
( 12002bl ); a forthcoming paper (Belczynski et al., in preparation) will discuss 
these changes and the evolutionary physics underlying them in significantly 
greater detail. 
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Like any population synthesis code, it evolves randomly chosen binaries from 
their birth to the present, tracking the stellar and binary parameters. For 
any class of events that is identifiable within the code, such as supernovae 
or DCO mergers, we estimate event rates by taking the average event rate 
within the simulation (i.e., by dividing the total number of events seen within 
some simulation by the duration of that simulation) and renormalizing by a 
scale factor that depends on properties of the simulation (i.e., the number of 
binaries simulated and the binary birth mass distributions assumed) and the 
Milky Way as a wh ole (i.e., the present-day star formation rate); see Eq. (2) 



and Appendix A of lO'Shaughnessy et al.l (l2005al ) for details. Specifically, our 



simulations^ are normalized to be consistent with an assumed Milky Way star 
formation rate of 3.5M yr _1 . 

When constructing our archived population synthesis results, we did not choose 
to record detailed information about the nature and amount of any mass 
transfer onto the first-born NS. We therefore cannot reconstruct precise popu- 
lation synthesis predictions for the NSNS(vw) formation rate. However, we do 
record whether some mass transfer occurs, and the nature of the mass transfer 
mechanism. The conventional mass transfer mechanism - dynamically stable 
Roche-lobe overflow - inevitably produces a disk around the compact object 
and can potentially spin that object up. Other mechanisms, such as (possibly 
hypercritical) common-envelope evolution, presumably involve a substantially 
more spherical accretion flow; the specific angular momentum accreted may 
be substantially lower, possibly not enough to recycle the neutron star. Thus 
for the purposes of identifying a class of potentially recycled ( "visible" ) wide 
NS-NS binaries, we assume any system which underwent mass transfer (dy- 
namically stable in the case of wide systems) produced a recycled NS primary. 



5.4-2 Fitting results 



As in previous studies such as lO'Shaughnessy et al.l (j2005al ). here we have cho- 



sen to vary seven (7) model parameters in the synthesis calculations. These 
choices are strongly guided by our past experience with double-compact-object 
population synthesis and represent the model parameters for which strong de- 
pendence has been confirmed. These seven parameters enter into every aspect 
of our population synthesis model. One parameter, a power law index r in 



Our approach gives only the average event rate. The present-day merger rate 
agrees with this quantity when most mergers occur relatively promptly (e.g., < 100 
Myr) after their birth. Some DCOs - notably double BH binaries - have substan- 
tial delays between birth and merger, introducing a strong time dependence to 
the merger rate. The technique described above will significantly underestimate 
these rates. This point will be addressed in considerably more detail, both for the 
Milky Way and for a hetero geneous galaxy population in a forthcoming paper by 



O'Shaughnessv et al.1 (|2006al h 
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Fig. 11. Probability distributions for merging (top) and wide (bottom) NS-NS for- 
mation rates, as predicted by population synthesis. The light shaded region shows 
the 95% confidence interval consistent with observations of binary pulsars, as de- 
scribed in Sec. 5.1. Because of sparse sampling and relatively poor data, our fit for 
the formation rate of visible, wide NS-NS binaries is significantly more uncertain 
than other fits. 

our parameterization of the companion mass distribution, influences the ini- 
tial binary parameter distributions (through the companion masses). Another 
parameter w, the massive stellar wind strength, controls how rapidly the mas- 
sive progenitors of compact objects lose mass; this parameter strongly affects 
the final compact-object mass distribution. Three parameters V\ t 2 and s are 
used to parameterize the supernovae kick distribution as a superposition of 
two independent maxwellians. These kicks provide critical opportunities to 
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Fig. 12. Probability distributions for merging (top) and eccentric (bottom) WD-NS 
formation rates, as predicted by population synthesis. The light shaded region shows 
the 95% confidence interval consistent with observations of binary pulsars, as de- 
scribed in Sec. 5.1. The dark shaded region ex tends this constraint i nterval by an 
estimate of the systematic error in each fit; see lO'Shaughnessv et al.1 (j200fiah . 



push distant stars into tight orbits and also to disrupt potential double neu- 
tron star progenitors. Finally, two parameters a\ and f a govern energy and 
mass transfer du r ing cer tain types of binary interactions; see Section 2.2.4 of 
Belczynski et al.l ( 12002bl ) for details. To allow for an extremely broad range 
of possible models, we conside r the specific parameter ranges quoted in §2 
of lO'Shaughnessv et al.l (l2005al ) for all parameters except kicks; for supernova 
kick parameters, we allow the dispersion t> lj2 of either component of a bimodal 
Maxwellian to run from to 1000 km/s. 
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Fig. 13. Probability distributions for SN Ib/c (top) and SN II (bottom), as pre- 
dicted by population syntheses with StarTrack. The light shaded region shows 
the (~ 2a) interval con sistent with observational constraints of supernovae, from 
Cappellaro et al.1 (jl999T ). A dark shaded region (not visible in this plot) indicates 
the se constraints, augmented by an estimate of the systematic errors in our fits; 
O'Shaughnessy et al.l ( 2006al ) . These constraints provide no information about 



sec 



population synthesis. 

Since population synthesis calculations involve considerable computational ex- 
pense, in practice we estimate the merger rate we expect for a given combi- 
nation of population synthesis model parameters via seven-dimensional fits to 
an archi ve of roughly 3000 detailed sim ulations, as presented and analyzed in 
detail in O'Shaughnessy et al. ( 2006al ) (hereafter OKB). However, these fits 
introduce systematic errors, which have the potential to significantly change 
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the predicted set of constraint-satisfying models. For this reason, OKB also 
estimate the rms error associated with each seven-dimensional fit. Finally, 
OKB demonstrate that, by broadening the constraint-satisfying interval by 
the fit's rms error, the predicted set of constraint-satisfying models includes 
most models which actually satisfy the constraints. For this reason, the dark 
shaded regions in Figures 11, 12, and 13 account for both observational un- 
certainty in the Milky Way merger rate and for systematic errors in the fits 
against which these observations will be compared. 

5.5 Prior distributions 

Lacking knowledge about which population synthesis model is correct, we as- 
sume all population synthesis model parameters in our range are equally likely. 
A monte carlo over the seven-dimensional parameter space allows us then to 
estimate the relative likelihood, all things being equal, of various DCO merger 
rates (shown in Figures 11,12) and supernova rates (Figure 13). Also shown 
in these figures are the observational constraints described in Sections 5.2.1, 
5.2.2 and 5.3 (shown in shaded gray). Finally, the dotted lines in Figure 15 
show the a priori population synthesis predictions for BH-BH, BH-NS, and 
NS-NS merger rates in the Milky Way. 

5.6 Applying and employing constraints 

Physically consistent models must reproduce all predictions. Supernovae rates, 
being ill-determined due to the large observational systematic errors men- 
tioned in Sec. 5.3, are easily reproduced by almost all models at the 2a level; 
see Fig. 13. However, population synthesis models do not always reproduce ob- 
served formation rates for DCOs, as shown in Figs. 12 and 11. These four 95% 
confidence intervals implicitly define a (0.95) 4 ~ 81% confidence interval in 
the seven-dimensional space, consisting of 9% of the original parameter volume 
when systematic errors in our fitting procedure are included. In other words, 
we are quite confident (80% chance) that the physically-appropriate parame- 
ters entering into StarTrack can be confined within a small seven-dimensional 
volume, in principle significantly reducing our model uncertainty. 

In Figure 14 we show one-dimensional projections onto each coordinate axis 
of the constraint-satisfying volume - in other words, the distribution of values 
each individual population synthesis parameter can take. For the purposes 
of this and subsequent plots, we use the following seven dimensionless pa- 
rameters Xk that run from to 1: X\ = a/3; £2 = w; X3 = v i/(1000km/s), 
X4 = i>2/(1000km/s); x$ = s; xq = a\, and x^ = f a . As seen in the top 
panel of this figure, the conditional distribution of supernovae kicks bears 
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a surprising resemblance to the observed p u lsar kick distribution (see, e.g. 
Arzoumanian et aD 19991 ; Hobbs et al. 2005 ; Faucher-Giguere fc Kaspi 2006 



and references therein) , even though no information about pulsar motions have 
been included among our constraints and priors. The population synthesis mass 
transfer parameter f a is also well-constrained, with a strong maximum near 
x — 0, implying that mass-loss fractions of about 90% or higher are favored in 
non- conservative, but stable, mass transfer episodes. Also, common envelope 
efficiencies aX of about 0.2-0.5 appear to be favored by the constraints. The 
rest of the one-dimensional parameter distributions are nearly constant and 
thus uninformative. Though small, the constraint-satisfying volume extends 
throughout the seven- dimensional parameter space. 



Evaluating fits for other DCO merger rates on the constraint-satisfying model 
parameters provides revised estimates for various phenomena of interest, such 
as for the BH-BH, BH-NS, and NS-NS merger rates. Figure 15 shows smoothed 
histograms of the BH-BH, BH-NS, and NS-NS merger rates, both before 
and after observational constraints were applied. Additionally, a smoothed 
curve shows the results modulated by the effect of systematic error. These re- 
sults sho uld be contrasted with the dis tinctly different post-constraint results 
shown in lO'Shaughnessy et al. fl2005bh (cf. their Figure 5 and C7). Notably, 
even though two fewer constraints were imposed, they find a significantly 
smaller (2% versus 9%) constraint-satisfying volume. In l arge part, these differ- 



ences c an be ascribed to including systematic errors: in lO'Shaughnessy et al 



(J2005bJ), preliminary data for the sparse and poor-quality NSNS(vw) sample 
was fit and applied without any account for fit-induced errors. However, our 
calculation also differs at several fundamental levels from the original approach 
: (T) the space of models stud ied is larger, covering more area in v%, v 2 (see, e.g., 
O'Shaughnessy et al.ll2006ah : (ii) the observed sample of merging NS-NS bi- 
naries is compared with the total merging NS-NS formation rate predicted 
from population synthesis, rather than with subset of merging NS-NS which 
undergo mass transfer; and critically (iii) the constraints are more numerous. 



Finally, in Figure 16 we estimate the range of initial and advanced LIGO de- 
tection rates implied by these constrained results. This estimate combines 
our constrained population synthesis results; our estimates for systematic 
fitting error; the approximate LIGO range as a function of chirp mass M c : 
1 ! 6 where d Q 



d = d (A4 C /1.2M Q ) where d = 15y3/2, 300 Mpc for initial and advanced 
LIGO, respectively; estimates for the average volume detected for a given 
species, as represented by ( m} 5 / 6 ) = 111, 5.8, 2M^ /6 for BH-BH, BH-NS, 
and NS-NS, respectively (see lO'Shaughnessy et al.l l2005bl ); and a homoge- 



neous, Euclidean universe populated with Milky Way-equivalent galaxies with 
density O.OlMpc -3 , each forming stars at rate 3.5M yr _1 . Updated rates that 
address a number of these simplificat ions will be examined in tw o forthc oming 
papers, lO'Shaughnessy et al.l (]2006a ) and lO'Shaughnessy et al.l (]2006bl ). 
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Fig. 14. Differential probability distributions dPk/dx denned so Pk(x) is the fraction 
of all models consistent with all DCO observations and with Xk < X. The top 
panel shows the distributions for the 3 kick-related parameters X3,X4,xs (dashed, 
solid, and dotted, respectively, corresponding to vi, t>2, and s). The bottom panel 
shows the distributions for x\ (the mass-ratio distribution parameter r, appearing as 
the solid nearly constant curve) , and the three binary-evolution-related parameters 
X2,X6,X7 (dashed, dotted, and solid, respectively, corresponding to w, aX, and f a ). 
The distribution of f a = x-j exhibits a strong peak somewhere between — 0.1; our 
choice of smoothing method causes all projected distributions to appear to drop to 
zero at the boundaries, as it involves averaging with empty cells. 
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Fig. 15. A plot of the probability distributions for the BH-BH (top), BH-NS (center), 
and NS-NS (bottom) merger rates per Milky Way-equivalent galaxy, allowing for 
systematic errors in the BH-BH, BH-NS, and NS-NS fits. The solid curves result 
from smoothing a histogram of results from a random sample of population synthesis 
calculations with systematic errors included. The dashed curve shows the same 
results, assuming model parameters are restricted to those which satisfy all four 
DC O constraints (both WD-NS a nd both NS-NS constraints). Compare with Figure 
5 of lO'Shaughnessv et al.1 (j2005bh . 
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Fig. 16. Range of expected detection rates for initial (top) and advanced (bottom) 
LIGO for BH-BH (dashed), BH-NS (solid), and NS-NS (dotted) binaries, based on 
m odels that satisfy obs e rvation al constraints and under the assumptions outlined 
in Q'Shaughnessv et al. ( 2005bl ) (e.g., using a fixed chirp- mass spectrum and LIGO 
range for detection, as well as assuming Milky- Way like galaxies fill the universe 
with density O.OlMpc -3 ). These estimates incorporate systematic errors due to the 
fitting procedure. 

5. 7 Conclusions 



Q'Shaughnessv et al. fl2005bl ) performed a proof-of-concept calculation to com- 
pare population synthesis simulations with observations, interpret the re- 
sulting constrained volume, and apply the constraint-satisfying parameters 
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to revise a priori predictions of astrophysically critical rate results. Here we 
demonstrate that even when systematic model fitting errors are aggressively 
overcompensated for, observational constraints still provide information about 
population synthesis parameters. These results include both surprisingly good 
agreement with pulsar kick distributions and strong constraints on at least one 
parameter involved in binary evolution. Since in almost all cases our system- 
atic errors appear much smaller than those from observations, we are confident 
that now observational limitations, rather than computational ones, primarily 
limit our ability to constrain population synthesis results. In particular, we ex- 
pect stringent tests of population synthesis models are possible, beginning by 
imposing more observational constraints than varied model parameters (seven 
in our case). 

Our approach also strongly suggests that population synthesis estimates for 
DCO merger rates vary in a very complex manner, as indicated by the large \ 2 
produced by even th e best low-order polynomial fits to often often very high- 



quality data (see, e.g. lO'Shaughnessy et al.ll2006al ). We strongly suspect that a 



more physically motivated class of basis functions could achieve a much better 
fit; we intend to explore this possibility, along with the reliability of purely 
nonparametric fits, in a future study. 

Our analysis remains predicated upon a correct identification of all parameters 
and input physics critical for the formation of double compact objects. Ad- 
mittedly, our understanding of binary evolution continues to evolve; however 
multiple population studies from different groups over the years lead to very 
similar conclusions about the basic input physics that primarily affects the for- 



mation rates of double compact objects (IBelczvnski et al.ll2002bl ; iHurley et al 



20021 ; iPortegies Zwart fc Verbuntlll996l ; iFryer et al.lll998l ) . In this paper we use 



a relatively new vers i on of t he StarTrack code that is described in great detail 
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Belczynski et al.l (l2006bl ). A forthcoming paper (Belczynski, Kalogera, & 
Bulik, in preparation) will discuss how these updates affect the evolutionary 
history and formation rates in significantly greater detail. 

Our analysis also remains predicated upon a correct and complete interpre- 
tation of the observational sample and associated systematic and model un- 
certainties. In some respects, however, our models for observations and their 
biases may be incomplete or not representative. We effectively assume pulsar 
recycling is required to detect binary pulsars in our implicit hypothesis that 
only one pulsar in NS-NS systems, the one observed, was ever likely to have 
been detected. And finally we use a canonical value for pulsar beaming cor- 
rection factor f a for several systems with known spins but unknown beaming 
geometry. (We hope to address these questions in a forthcoming paper.) 

Last, our analysis of DCO merger rates remains predicated on the assumption 
that the Milky Way formed stars at a relatively uniform rate for the last 
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10 Gyr. However, as will be demonstr ated at greater length in two forthcoming 



papers (lO'Shaughnessy et al.l l2006bl ldj. if the star formation rate increased 



significantly when the Milky Way was young, then since most binaries would 
then form early in the lifetime of the Milky Way, the young-galaxy contribution 
could conceivably overwhelm the present-day contribution. In other words, the 
present-day Milky Way merger rate for BH-BH binaries, for example, could 
be completely determined by the star formation rate nearly 10 Gyr ago. 

And to be comprehensive regarding the factors neglected here, our population 
synthesis parameter study has not varied over all p lausible models. We have 



neither considered purely polar supernova kicks (see iJohnston et al.ll2005l ) nor 
changed the maximum NS mass from 2M nor even allowed for a heteroge- 
neous birth population of single stars and binaries (a binary fraction of 100% 
is assumed in all models). And of course we continue to use a single fixed 
star formation rate for the Milky Way, rather than treat it as an independent 
uncertain but constrained parameter. 

However, all these limitations can be remedied by closer study. Thus, we ex- 
pect that either (i) all DCO merger rates will become determined to within 
0(50%) or better (based on the comparative accuracy to which we know NS- 
NS merger rates), or (ii) experimental data will conflict with the prevailing 
notion of binary population synthesis, revealing flaws and limitations in clas- 
sical parametric models for binary stellar evolution. 



6 Open Issues 



It has now been more than 30 years since the discovery of the first binary 
pulsar that happened to also be the first binary with two neutron stars. This 
discovery has essentially led to the development of a whole sub-field in compact 
object astrophysics and progress in this area has also greatly contributed to 
the progress in the development of gravitational- wave observatories over these 
years. 

As we anticipate the direct detection of gravitational waves from binary inspi- 
rals and the possible, associated discovery of black hole binaries, we also take 
advantage of pulsar timing observations to the fullest and continually try to 
improve our understanding of the origin of tight double compact objects. It 
is clear that a lot of progress has been particularly made in the last few years 
with the discovery of the highly relativistic, first double pulsar and additional 
double neutron stars, and the tremendous progress in proper motion and pulse 
profile evolution measurements. More interestingly, our standard way of think- 
ing about double compact object formation is called to question with a number 
of intriguing suggestions and counter-suggestions put forward: do neutron stars 
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form through two different physical mechanisms dep ending on the mass of their 



immediate progenitor at the time of core collapse (IPodsiadlowski et al.l 12005 



Piran & Shaviv 


2005a: 


Willems et al. 


2006; 



majority of neutron stars receive minimal supernova kicks with magnitudes sig- 
nificantly smaller than SOkms -1 ( van den Heuvellbooi IChaurasia and Bailed 



20051 : llhm et al.ll2006l : iDewi et al. 



20051 ) . Furthermore, some of the older ques- 



tions still persist: what is the basic physical mechanism for core collapse that 
also drives supernovae? Do most black holes form through fallback onto proto- 
neutron-stars and supernovae or through direct collapse? What is the mass 
function of black holes and how does it depend on metallicity? Why are neu- 
tron star masses so close to one another in double neutron stars? What really 
determines the birth spins and spin evolution of neutron stars and black holes? 

Starting with the very basic question of the core collapse physical mechanism 
Hans Bethe was intrigued by many of these questions and naturally his inter- 
est evolved towards binary compact objects, their evolutionary history, and 
their potential as gravitational wave sources. Given the progress we have expe- 
rienced since the discovery of the first DNS binary, it is reasonable to expect 
that with the combination of developments in pulsar searches and observa- 
tions, of the rise of gravitational-wave detections and astrophysics, and of the 
progress in theoretical physical understanding and computational astrophysics 
we will be able to answer some of the basic remaining questions raised by the 
existence of binary compact objects in nature. 
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